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Abstract 

Front form dynamics is not a manifestly rotational invariant formalism. 
In particular, the requirement of an invariance under rotations around 
the transverse axes is difficult to fulfill since the corresponding opera- 
tors are complicated and involve the interaction. All approximations 
in solving Quantum Field Theories using front form dynamics, as in- 
evitable as they are, are destined to destroy full Poincare invariance. 
In the present work it is investigated, to which extent rotational invari- 
ance is restored in the solution of a light-cone quantized field theory. 
The positronium spectrum in full (3-1-1) dimensions is calculated at an 
unphysically large coupling of a = 0.3 to be sensitive for terms break- 
ing rotational invariance and to accustomize the theory to future QCD 
applications. The numerical accuracy of the formalism is improved to 
allow for the calculation of mass eigenvalues for arbitrary components 
of the total angular momentum. We find numerically degenerate eigen- 
values as expected from rotationally invariant formalisms and the right 
multiplet structures up to large principal quantum numbers n. The re- 
sults indicate that rotational invariance is unproblematic even in front 
form dynamics. Another focus of the work relies on the inclusion of the 
annihilation channel. This enlargement of the model is non-trivial and a 
consistency check of the underlying theory of effective interactions. The 
correct numerical eigenvalue shifts, and especially the right hyperfine 
splitting are obtained. Moreover, the cutoff dependence of the eigenval- 
ues is improved drastically by the annihilation channel for the triplet 
states. The implications of the applied effective Hamiltonian approach 
are discussed in detail. 
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Introduction 



After a long line of efforts, Quantum Chromodynamics [0, H, ||, H emerged around 
1973 as the theory that correctly describes the interaction between quarks (fermions), 
considered to be elementary, through the action of the non-abelian gauge field of gluons 
(bosons). The hierarchy of particles seems to be the following. The hadrons, that is 
baryons (nucleons) and mesons, consist of quarks. The gluons produce the interaction 
between quarks. They provide the force, a strong analogue of the van-der-Waals force, 
that allows nucleons to bind together to form the atomic nucleus. This point of view 
concerning the nuclear forces clarified the former picture that the intermediate boson 
field is that of the mesons 0. 

An open question in modern theoretical physics is the link between the high and 
the low energy domains of Quantum Chromodynamics. The description of the high 
energy region, which was inspired by experiments in the late sixties, is well established. 
Contrary, the low energy region of the theory cannot be attacked by standard per- 
turbation theory because the coupling depends on the scale and grows large in this 
region. Instead, this region must be described by more phenomenological models such 
as the constituent quark model, or chiral perturbation theory. A major problem in this 
undertaking to connect the high and the low energy regimes is the running coupling 
constant of QCD which depends very strongly on the four momentum transfer between 
the interacting particles. Moreover, it grows large for small momentum transfer so that 
perturbation theory breaks down at a certain point. An important task of theoretical 
high and medium energy physics is therefore the construction of methods that can pro- 
duce hadronic spectra and wavefunctions. The importance of the hadronic properties 
and their wavefunctions for physical observables, such as structure functions and decay 
constants, is undeniable. Moreover, their calculation would allow for rigorous tests of 
assumptions necessary in models of nuclei and hadrons. 

The history of methods to calculate quantum mechanical bound states is as old 
as Quantum Mechanics itself 0. But in hadronic physics the large coupling implies 
that bound particles form highly relativistic systems. This fact enforces the rigorous 
application of covariant field theories. Non-relativistic potential models [|^, ^, ^ are 



only applicable for heavy systems, as for instance the J/ip |T^, |TT| or the T |T2 
Other covariant methods are lattice gauge theory or the evaluation of the Bethe- 
Salpeter equation |T^. Within lattice gauge theories only the ground state or the 
first few excited states are tractable and structure functions are extremely difficult 



to observe directly [|T5|. The application of the Bethe-Salpeter equation||T^ suffers the 
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typical fate of strongly bound systems: The wavefunctions are dominated by relativistic 
momenta and contain a considerable number of quark- ant iquark pairs. Therefore, even 
the vacuum is very complicated. 

The latter point is problematic for all quantum field theories, if a theory is quan- 
tized at equal usual time. In light-cone quantization, the results are frame indepen- 
dent, as opposed to equal time quantization. There, the Lorentz boosts contain the 
interaction and are therefore complicated dynamical operators, whereas in light-cone 
quantization all Lorentz boosts are kinematic [1^, Chap. 2.B]. Consequently, light-cone 
quantized field theories seem to be a promising tool to understand the physics of rel- 
ativistic QCD-bound states. However, neither have all problems (Poincare invariance, 
zero modes, etc.) of this Hamiltonian approach been solved, nor is it clear if and how 
effective theories [|18], |19|, ^ which have been proposed to deal with the problems asso- 
ciated to finite cutoffs, will work. Some of these deficiencies are due to the neglection 
of Hamiltonian field theory in favor of the action-oriented approaches of Feynman, 
ET. AL |2^, 23]. In any case, it is better to establish that a new method works for 



known problems before applying it to yet unsolved problems. Bearing this in mind, the 
basic motivation for the present work is twofold. On the one hand, there is the obvious 
obstacle that light-cone quantized theories are rotationally not manifestly invariant by 
construction. One main part of this thesis is consequently dedicated to the investiga- 
tion of the associated question, what this means to the results of such a theory: do the 
results obey rotational symmetry (as observed in nature) or not? Because we have to 
compare to other calculations and methods, we choose a realistic, well-understood, but 
non-trivial problem, namely the QED-bound state par excellence: positronium. Our 
choice has another cause. The formahsm to be applied is ideally suited for this case 
of two particles of same mass. In ordinary approaches, the severe problems of recoil 
corrections make it much easier to calculate the hydrogen or muonium spectra than to 
solve for the eigenvalues of positronium. 

On the other hand, there is the need for an effective formalism for gauge theories: 
nobody ever has solved rigorously a relativistic many-body theory in one time and 
three space dimensions (3+1 dimensional theory). It is therefore necessary to use such 
an effective approach and to investigate how reasonable its results are. We have done 
so. But after the successfully applying of the effective method, we go one step further 
in asking how the excellent results can be understood within a much more general 
formalism, originating in the Lagrangian density of a gauge theory. 



Outline of the present work 

This work is structured as follows. We start in Chapter ^ with an overview of the 
development and difficulties of the approach to field theories in front form dynamics. 
In Chapter ^, the positronium model used in the present work is introduced. Starting 
with a Hamiltonian field theory quantized on the light-cone, the steps are described that 
lead to the formulation of an effective integral equation in momentum space. Although 
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the derivation of tlie Hamiltonian involved the apphcation of the DLCQ formahsm, 
described in the following section, we stress the fact that we use only the continuum 
formulation of QED(3+i) in front form dynamics in the present work. Discretizations 
occur only in the context of numerical integrations, which have nothing to do with 
physical considerations. Errors due to numerical artifacts, such as effects of a finite 
number of integration points A^, are discussed. Cutoffs determined by physics emerge 
in two aspects. One is the problem of divergencies occuring because transverse mo- 
menta can become infinite. Therefore, one must introduce a regulator. We analyze the 
results of our model concerning their dependencies on this cutoff. The second is the 
restriction of sectors with different particle content, sometimes referred to as Tamm- 
Dancoff approach. We will show that (a) a theory truncated in this way is not solvable, 
(b) there is a possibility to treat the associated singularity by adjusting the form of the 
effective interaction, and (c) the treatment of the singularity can be understood within 
a general formalism of effective interactions. The latter point is discussed in Chapter |^. 

We improved the numerical methods which were used to solve the effective inter- 
action before [|^. This leads to an improved convergence of the eigenvalues with the 
number of integration points. However, the basic intention is the need of this improve- 
ment to attack the generalization of the model, described in Chapter ^ More effort 
than in is put in the investigation of higher excited states to show the range of 
applicability of the model. We find the correct multiplet structure up to a principal 
quantum number n=5. Detailed tables allow for the comparison of the results of the 
present work with those of equal-time perturbation theory. We emphasize the use of 
a huge coupling constant, a=0.3, to test our non-perturbative approach: it is approxi- 
mately 40 times larger than the physical coupling constant of QED. In the calculation 
of the hyperfine splitting, we point out that the comparison of our results with per- 
turbation theory is problematic, because results of the latter depend noticeably on the 
order of the calculations at large couplings. 

As a main point of the thesis, we generalize the formalism of Chapter ^ to arbitrary 
z-components of the total angular momentum, J^. This is done in Chapter ^ Before, 
calculations were performed in the Jz=0 sector of the theory only. To interpret the 
results, the theory of the Poincare group is studied in the introductory section of this 
chapter in the context of front form dynamics. We find that states with the same 
total angular momentum J, but different Jz, are numerically degenerate. This result is 
expected from physical considerations, but surprising because the associated operators 
of transverse rotations are dynamical in front form dynamics although they do commute 
with the light-cone Hamiltonian. As a consequence of our calculations, we can classify 
the positronium eigenstates according to their quantum numbers of and Jz, just as in 
the equal-time formalism. The wavefunctions are analyzed concerning their symmetry 
properties. Besides, we show that the wavefunctions are not rotationally invariant in 
any sector of Jz- 

After the application of effective interactions and the generalization of the theory 
to arbitrary Jz, the next crucial point is the introduction of the annihilation channel 
into the theory (Chapter We enlarge the Fock basis by the one-photon state. 
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typical for QED and absent in QCD. It is shown that the degeneracy of states is 
maintained. This is astonishing, because we find a total separation of the dynamic and 
instantaneous diagrams involved in the different Jz sectors. Another consequence of the 
implementation of the annihilation channel is the stabilizing of eigenvalues, concerning 
their dependence on the cutoff A. 

Chapter ^ deals with the problems constructing effective interactions in front form 
dynamics. We show that the prescription to treat the singularity attached to a trun- 
cation of the Fock space is a consequence of the formalism of effective interactions. A 
summary and a discussion of our results follow. 

All basics and technical details are consequently cast into the appendices to achieve 
a clear and structured line of arguments in the main part of the thesis. We address 
the attention of the reader especially to Appendix |B| where QED on the light-cone 
is described, and to Appendix |y where the numerical methods are explained. The 
latter should be studied when details of calculations seem unclear in the main text. In 
Appendix |^ the applied renormalization scheme is described. The other appendices 
mainly supply necessary formulae. The source code of the computer program developed 
in the present work is listed in Appendix 0. 

Methodological sketch 

This section describes in short the methods applied in the thesis and their subtle points. 

• The underlying picture of a composed system in front form dynamics is that of 
a state consisting of valence particles carrying its outer quantum numbers, and 
of arbitrarily many virtual particles, be it gauge particles or fermion-antifermion 
pairs. A state is written in Fock-space representation as 

1^) = ^jV^s) + ^otI^ot) + 1^^99-7) + '^mggl'^q'igg) + ■■■■ 

Here, g denotes a gauge boson (photon or gluon) and q a fermion (electron or 
quark). In QCD, the first term in this expansion is absent because it is not color 
neutral. 

• To construct the Hamiltonian, the following steps are performed. First, starting 
with the Lagrangian density, the generators of the Poincare group are calculated.^ 
Next, the dependent fields are expressed in terms of the dynamic fields. 

• The independent fields are expanded in plane waves and are discretized by impos- 
ing periodic or anti-periodic boundary conditions along all space-like directions 
in a given volume. 

^In actual calculations, one restricts oneself to calculating the momenta P^, because the angular 
momentum tensor M^'^ is irrelevant for the mass spectrum. Sometimes Jz is constructed. 
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It is important that the expansion into plane waves takes place in momentum 
space. By this, assuming a reasonable choice of boundary conditions, all points in 
this space are equivalent and the momentum of the center of mass is well-defined. 
The enormous improvement of convergence achieved by such a prescription is seen 
already in non-relativistic many body theory . 



When working in the light-cone gauge, ^4+ = 0, the theory is consistent in the 
normal mode sector only (cf. next section). The zero modes are important for the 
vacuum structure |2^, but most probably not for bound state spectra. 



Having quantized the theory by postulating canonical commutation relations for 
the creation and destruction operators, the Poincare generators are written as 
functionals of these operators. 

The reference state, "the vacuum", is an eigenstate of the Hamiltonian with 
eigenvalue zero. This follows from the fact, that creation of particles out of the 
vacuum is forbidden because of the conservation of longitudinal momentum. 

Spectrum and wavefunctions are obtained by solving the eigenvalue problem 

where the ipn are eigenfunctions to the mass (squared) eigenvalue M„. This 
matrix equation and its solution by diagonalization are the endpoint of the so- 
called DLCQ formalism. We discuss its difficulties in the next chapter. 

We can perform the continuum limit of this matrix equation. This is possible in a 
direct way, as opposed to the non-trivial continuum limit, for example, of lattice 
gauge theories. The matrix equation is mapped into an integral equation. 

The task is now to find an appropriate scheme to solve for the spectrum of this 
(infinite dimensional) operator. We use a restriction of the Fock space together 
with an effective interaction to account for the effects of the truncated states, and 
an explicit cutoff because of the unrestricted transverse momenta. To exploit and 
to justify this approach is the main part of the present work. 

The so-obtained effective integral equation is solved with the appropriate method 
of Gaussian quadratures, which leads in the end to the diagonalization of a finite 
dimensional matrix. 

The dependence of the solutions on the unphysical parameters (cutoffs, etc.) has 
to be investigated to test the significance of the obtained results. 
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Chapter 1 

Field theories in front form 
dynamics 



1.1 Historical survey 

Light-cone coordinates^ were introduced into the Hamiltonian field theory by DiRAC 
[P7| in 1949. According to his definition an operator is called a Hamiltonian if it 
propagates a physical system in a fixed direction in space-time. This direction is subject 
to certain constraints, but not unique. DiRAC listed three different possibilities^ in 
selecting such a "general time" . The most prominent ones are the usual time x° {instant 
form dynamics) and the light-cone time {front form dynamics) . This Hamiltonian 
formalism attracted little attention because of the spectacular achievements (e.g. [28, 
the anomalous magnetic moment of the electron, Eq. (1.112)]) of the action- oriented 
work of TOMONAGA [^, SCHWINGER and Feynman |2^, launched at the 



same time. 

The rejuvenation of the Hamiltonian method in field theory, in particular studying 
front form dynamics, came in several steps. First, Weinberg |2^, using (p^ theory in 
the infinite momentum frame, discovered that the creation and annihilation of particles 
out of the vacuum is forbidden for —>■ oo and that the corresponding divergencies 
are absent in this frame. After the derivation of the rules for front form perturbation 
theory |^ , it was an important success to show the equivalence of this theory with the 



Feynman rules of ordinary perturbation theorypl]][B3]. The treatment of non-abelian 



gauge theories and the description of exclusive QCD processes are described in |34 

To render the theory tractable and to finally implement it on a computer, the 
ambitious program of Discretized Light-Cone Quantization was proposed by Pauli 
and Brodsky in 1985. The formalism was applied first to the Yukawa model in (1 + 1) 
dimensions, then to scalar QED(i+i) followed by calculations of the theory 



33, and QCD(i+i) |^ Commencing with a paper by Eller, Pauli and Brodsky 



^Thc coordinate vector is a: = (a:+, x , xj_) with = x^ ± x^, xj_ — {xi,X2)- Cf. Appendix 
^Actually, there are five such quantization surfaces 
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P0| , there is a large literature discussing the (massive) Schwinger model |4^, |43|] . 

All of the above examples have been studied in 1+1 dimensions. Quantum Elec- 
trodynamics was the first realistic gauge theory in physical space-time to be treated 
with DLCQ H m, 



15| , where results were compared with experimental data or 
other theoretical work. Subsequently, an attempt was made to apply this method to 
QCD(3+i) Most recently, this was done to a large extent in the so-called collinear 
i.e. effectively in lower dimensions. Here, the method of the transverse 



model |8 



lattice 1 49, pOi seems to open a promising way to proceed, but is still under construction. 



Nevertheless, QED(3+i) is a "milestone" for non-perturbative approaches in light-cone 
quantization on the way to understand full QCD. QED contains most of the fundamen- 
tal difficulties found in QCD and can be attacked likewise in a truncated Fock basis. 



1.2 Problems of front form field theories 

The disadvantages using front form dynamics include a certain non-conformity of the 
inertial system used, especially the counter-intuitive missing of a well-defined angular 
momentum. As mentioned in the introduction, the operators of rotations around the 
transverse axes are complicated, i.e. contain the interaction. Consequently, although 
the rotation operator jTs is kinematic, the states of a system cannot be classified with 
respect to total angular momentum J^. Of course, this is not necessary from the out- 
set, and calculations have been performed without using equal-time quantum numbers, 
for example by restricting to a sector of a definite z-component of the total angular 
momentum J^. However, the physical results of a calculation should be independent of 
the mathematical method applied. A priori it is not clear in front form dynamics, if 
the results of a calculation will be rotationally invariant. It seems to be impossible to 
show analytically that these solutions do indeed obey rotational symmetry. This would 
require the diagonalization of an operator at least as complicated as the light-cone 
Hamiltonian itself. The problem, though fundamental for all Hamiltonian formulations 
of field theories, was not of primary interest before, since calculations were often per- 



formed using lower (typically 1+1) dimensional models. It has been shown [^, ^ that 
Lorentz covariance, and in particular rotational symmetry, is explicitly violated, if one 
evaluates front form perturbation theory at the one- and two-loop levels. Non-covariant 
counterterms have to be constructed to restore covariance. 

Because of its long term use, many phenomena have been investigated within equal- 
time formulation. The number of calculations is understandably much smaller in light- 
cone quantization. Much effort has been made therefore, to reproduce the results 
of instant form dynamics on the light-cone. A considerable difficulty in light-cone 
quantization is the role of the so-called zero modes. The zero mode of a function 



f{x , x±) in a fixed space direction y with interval length Ly is defined as (cf. e.g. 
Eq. (4)]) 

1 f^y 

{f{x))o := ^ / dyf{x,y), 
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where it is not integrated over the remaining space directions x. As a first step, the 
zero modes were omitted in actual calculations, because they are a set of measure zero 
among the modes of all fields. Their infiuence on the spectra of bound systems was 
considered to be negligible. 

However, it was found that the convenient light-cone gauge = is inconsistent 
with the inclusion of the zero modes into the formalism. In non-abelian gauge theories, 
for instance, the zero mode of the gauge field cannot be gauged away [^. It was 
often seen that by working in the light-cone gauge in the front form dynamics, spurious 



operators which are singular in the limit A;"'"=0 are created The renormalization of 



these ultraviolet (UV) and infrared (IR) divergencies was investigated in perturbation 



theory by Mustaki et al. [Q. Although the most divergent contributions are can- 
celled in the graphs of lowest order, the IR singularity remains. For example, the photon 
propagator picks up an additional l/fc"*" singularity. After regularization, it gives rise 
to difficulties only in higher orders, because then one must integrate over the photon 
momentum. A naive regularization via principal value runs into difficulties which can 
be overcome by the Leibbrandt-Mandelstam prescription [^. This method for 
Hamiltonian non-abelian gauge theories was derived by Basetto et al. ||56|| . 

To put the former results of calculations in light-cone quantization on a formally 
correct basis, the inffuence of the zero modes was examined systematically after prag- 



matically ignoring them. To name a few, Pauli, Pinsky and Kalloniatis E^, 



Werner, Heinzl et al. [^, and McCartor @ have made important contribu- 



tions to this subject. The problem of zero modes in QED is discussed by Kalloniatis 
and Robertson in [^. Recently, these more formal examinations have been accom- 
panied by numerical calculations. Vollinger |43| found in his investigation of the 
massive Schwinger model a vanishing inffuence of the zero modes concerning the spec- 
trum of this model. However, he considered a vacuum wavefunction that is independent 
of the variable 9, parameterizing the vacuum states. In literature the vacuum is 



found to be dependent on this parameter ( "^-vacuum" ) which can have an effect on the 
spectrum. As an important result, van de Sande |^ was able to show that DLCQ 
does give the correct linear increase of the squared mass eigenvalues with the fermion 
mass in this model. The work of Elser displayed a quadratic rise of this squared 
mass with the fermion mass. This effect was falsely attributed to the zero modes. The 
results of Elser were numerically unsatisfactory and the latter interpretation could be 
excluded in Ref |^^, where the method of van de Sande was used. It was proven 



that in the continuum limit the right slope is reached and is hidden by an extremely 
slow convergence typical of DLCQ. The solution to this problem resides in the use of 
counterterms. 

In conclusion, one has strong evidence for a small effect of the zero modes on the 
physics important for the spectra of bound systems. An application of the light-cone 
gauge, or equivalently, the consideration of the normal mode sector only, seems to be 
justiffed, at least in the QED case. Of course, in QCD there are large effects coming 
from the vacuum structure, such as chiral symmetry breaking. The origin of these 
effects has to be investigated separately. 
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1.3 Preceding work 

The investigation of QED within the DLCQ formahsm dates back to roughly 1988. In 
1990, Tang set up a matrix equation for QED(3_|_i) by truncating the Fock space to 
the sectors |ee) and |ee7). In order to solve the associated eigenvalue problem, he used 
the diagonalization of the discretized Hamiltonian, and variational methods. With help 
of the latter, he was only able to produce an upper limit for the mass-squared eigenvalue 
of the triplet ground state of positronium, three percent larger than the Bohr value at 
a = 0.6. But compared to the standard results of perturbation theory [Q, this upper 
limit is already below the well-known value for this triplet term. To produce significant 
results by diagonalizing the Hamiltonian matrix. Tang Ref. 1, p. 46], following his 
own estimates, would have had to include roughly 11 million Fock states for the large 
coupling constant a = 0.3. The main cause for this negative conclusion was the slow 
numerical convergence of the method applied. 

With the solution of the Coulomb problem in momentum space [pi], the method 



of Coulomb counterterms was introduced to improve numerical convergence. Kraut- 
GARTNER |Q applied this method to QED (3+1). He used an effective interaction, 
obtained from a projection of the |ee7)-sector onto the |ee)-sector. The correspond- 
ing effective integral equation was solved using Gauss-Legendre quadrature. His results 
show excellent convergence and coincide to a high degree of accuracy with the expected 
values. It is worth mentioning that, for the first time, not only the complete Bohr spec- 
trum is obtained, but also relativistic effects, like the hyperfine splitting, are correctly 
described within a non-perturbative approach in front form dynamics. 

Kaluza |45] applied the counterterm technique to the calculation of the spectrum 



of the light-cone Schrodinger equation. He calculated an eigenvalue for the positronium 
ground state that coincides approximately with the Bohr value at a = 0.3. For the 
calculation of relativistic spectra, Kaluza improved the diagonalization technique to 
enlarge the Fock space feasible with the computer equipment used. The convergence 
of his spectra is rather poor, because he did not use counterterms for the singularity 
of the relativistic problem. For a comparison with the light-cone Schrodinger equation, 
only the ground state is considered, which is certainly too large. 

WOLZ 1^ subsequently tried to solve the analogous problem for QCD. To include 



typical QCD effects (self-coupling of the gluons) he enlarged the Fock basis implement- 
ing the |gg(7(y')-sector. In order to solve for the eigenvalues, he applied a hybrid method 
by projecting the new Fock sector unto the other two sectors, considering the so-derived 
eigenvalue problem as a matrix equation. His results |H6| converge too slowly and suffer 



strong fluctuations as functions of the changing Fock space size. This is due to the fact 
that the counterterm technique included in the computer code is not apphed in the 
final calculations because the necessary two dimensional numerical integration are too 
time consuming. 
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Chapter 2 

The positronium model 



When quantizing a field theory on the hght-cone, it is convenient to work in the usual 
Fock basis. Consequently, one thinks of a composite physical system, e.g. a meson or 
a baryon, as consisting of fermionic valence particles, bosonic (virtual) gauge particles, 
and virtual particle-antiparticle pairs. Considering positronium as the physical system 
of special interest for the present work, the wavefunction reads: 

\pOSitronium) = + 1pee\lpee) + 1pee-f\lpee-y) + i'ee-y-,\'ipee^j) + ■ ■ ■ ■ 

The exact wavefunction is an infinite series and one has to impose simplifiying restric- 
tions in actual calculations. The obtained results have to be investigated concerning 
their dependencies on these restrictions. 

The Hamiltonian operators ruling the dynamics of a system can be derived using 
light-cone quantization. A rough sketch of the procedure is given in the introduction and 
its application to QED(3+i) is described in more detail in Appendix ^. It is important 
that in this work the continuum formulation of a light-cone quantized field theory 
is used. All discretizations come from numerical integrations, cf. Appendix 0. The 
operator 

is commonly called the invariant mass (squared) operator. For convenience I will refer 
to it as the light-cone Hamiltonian, although in the sense of DiRAC only P~ is a 



Hamiltonian. The light-cone Hamiltonian is obviously a Lorentz scalar, and so are its 
eigenvalues, which have the dimension of a mass squared. 

Our main task is to set up an effective, relativistic Hamiltonian operator, tractable 
either analytically or with the help of a computer, and to solve the eigenvalue equation 

H^cm=M'M, (2.1) 

whose solutions yield the mass (squared) eigenvalues and the corresponding eigenfunc- 
tions of our positronium model. The relevant matrix elements of this Hamiltonian are 
tabulated in Appendix y. The full Hamiltonian matrix elements for QED and QCD 
are given in MM and E^, , respectively. 
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2.1 The model 

To construct a manageable physical model for positronium, we proceed as follows. One 
is free to divide the Fock space into two arbitrary subspaces, called the P-space and 
the Q- space. Restricting the full Fock space to two sectors, being explicit, one with 
an electron positron pair and one with an additional photon, the associated projection 
operators onto these subspaces are defined to be 

P--=Y: l(ee)n)((ee)„| (2.2) 

n 

allQN 

and 

Q:= E I(ee7)„)((ee7)„|. (2.3) 

n 

allQN 

Within this limited Fock space, the model positronium cannot decay into photons, 
contrary to observation. Moreover, a virtual photon as an intermediate state in an 
interaction is impossible. The part of the hyperfine splitting connected to this anni- 
hilation graph is missing. The inclusion of this graph into the theory is the topic of 
Chapter ^ Nonetheless, the vector space is well-defined mathematically, and we can 
proceed to solve the eigenvalue equation ( |2.1|) . In our restricted Fock-space we have a 
2x2 block matrix 

= ( f ) . (2.4) 

\ -Hqp -Hqq J 

Of course, this truncation of the Fock space violates gauge invariance. It will be shown 
in Chapter ^ that one can treat the consequences of this violation within a general 
theory of effective interactions. 

The attempts to solve for the spectrum of this operator in a matrix equation P^, Q 



were bound to fail, as described in Chapter H because of computer capacity limitations. 
We apply a projection method (also used in many body theory^) to eliminate the Q- 
space at the expense of an in general more complicated effective Hamiltonian operating 
in P-space only. Formally, one can describe the projection as the solution of a system 
of coupled linear equations. The last of the two equations of the eigenvalue problem 



Eq. (|2.1| ) with the Hamiltonian (|2.4| ) reads explicitly 

HQp\ee) + HQQ\ee-i) = M^\en). 

We therefore solve for the state Q\ip) = {eej) and express it with the help of an inverse 
Hamiltonian, or resolvent: 

G{io) ■.= Q{uj-Hi^c)-'Q, (2.5) 

and obtain 



^It is described as mathematical method in [ |66| and appHed to problems in nuclear theory by Tamm 
IgtII and Dancoff isl. 
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1-x,-k^ 1-x',-k/ 

Figure 2.1: The seagull graph in P-space. 



x,k^ x',k'^ 




1-x,-kx 1-x'.-k' 



Figure 2.2: The iterated vertex interaction 
(x > x') in P-space. 



The a priori unknown eigenvalue was substituted by uo, referred to in the remainder 
of this work as redundant parameter. This parameter has to be fixed by an additional 
condition. 

Which interactions can occur in the effective |ee)-sector? The only possible graph 
in P-space is shown in Figure Because of the projection, one also has iterated 

graphs, like the one displayed in Figure (P^). The Q-space seems quite complicated 



at first glance, but only one graph, the lower right one of Figure ( |2.3| ), survives the 
so-called gauge principle of DLCQ formulated by Tang et al. [0, Ref. 2]. The main 
idea is to restore gauge invariance at tree level, which was destroyed by the Fock space 
truncation. The recipe is to omit all graphs which have intermediate states that are 
not contained in the (restricted) Fock space. Instantaneous particles are counted as 
real particles in this procedure. 



Still not all non-diagonal graphs in Q-space vanish. To simplify the problem, in 
the remaining seagull graph was omitted on an ad hoc basis. We shall show in Chapter 
^, using the formalism set up in [^], that this omission is not an assumption, but is a 
natural consequence of the projection mechanism. We arrive at the non/mear equation 

iJ£^M|^„M)=M2(.;)|^„H), 
where the states, the effective Hamiltonian 

Hfc{uj) := PH^cP + PHi^cQiuJ - H^c)-'QH^cP (2.6) 

and the mass eigenvalue depend on the redundant parameter u. In principle, we 
also have to satisfy the obvious constraint that 

Although the fixing of the parameter u is explained in Chapter ^, we introduce a 
method of calculating an analytic expression for this parameter. This method is the 
way the so-called uj* -trick was introduced before the work of Pauli on the Method of 



Iterated Resolvents [20|. The derivation given here should be considered as a plausibility 
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argument rather than as a strict proof. Nevertheless, it gives some intuition concerning 
the physics behind the procedure. 

A word on notation seems in ordeiQ. We label the fermion mass with ruf, the 
longitudinal momentum fraction with x, the transverse momentum with k±, and the 
helicity of a particle with A. The corresponding quantum numbers after an interaction 
have a prime (x', k'^] A'). 

We have seen that the energy denominator contains the unknown (parameterized) 
eigenvalue uj of the whole eigenvalue problem and the sector-Hamiltonian operating in 
Q-space. The latter consists of a kinetic part Mg, i.e. the free mass (squared) of the 
Q-space, and an interaction Vq: 

Hlcq = Ml + Vq. 

As pointed out above, this poses a difficult problem: a non-diagonal operator has to be 
inverted and as a constraint it has to be guaranteed that the masses should be equal to 
the mass parameter. The latter results in a mathematical fixpoint equation. To avoid 
at least the first of these difficulties, one can divide the interaction into a diagonal part 
(Vq) and a non-diagonal part 5Vq 

Vq = {Vq) + 5Vq. 

By defining formally 

T*:=uj- {Vq) = cI, c e R, (2.7) 

I changed the notation of |^ Eq. (5.9)] to stress two points. Firstly, T* is not just a 
fixed value of u. This would be inconsistent, because o; is a real number, whereas T* 
is a function of the light-cone momenta. Secondly, T* is a kinetic energy according to 
Eq. ( |2.7| ), where a potential energy is subtracted from a total energy. One can expand 
the resolvent around the diagonal interaction (Vq) 

(2.8) 



u-Ho T*-M?,- 6V 



Q 



T* -Ml T* -Ml ^T* -Ml- 5Vq ' 

As a first approximation we consider the first term of the expansion only. It does not 
contain any non-diagonal terms and is a simple c-number. We label 

V{x,x'-T*) ■=\x-x'\{T* - Ml). (2.9) 

This approximation has severe consequences: a coUinear singularity occurs. It is 
proportional to 

1 A(x,fcx,x',A;^;T*) 
V \x — x'\ 

^Cf. Appx.0. 
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-fV\/WV\/\ 



■VN/WVX/vi- 




Figure 2.3: Graphs of the instantaneous interaction in Q-space. 



We identify now the |ee)-sector with the P-space and the |ee7)-sector with the Q-space 
to calculate the function A by evaluating the two graphs Fig. (|2.1| ) and Fig. (|2.2|) 
according to the rules of front form perturbation theory given in Refs. 133, BH], but 



using the expression in Eq. (|2 
obtains 



for the summation over the intermediate states. One 



667 



UJ 



X 



X' 



where := (/Cg — ke)^ and Ig := {kg — k'^)^ are the momentum transfers of the electron 
and positron. One can use this singularity induced by the truncation of the series (|2.8|) 
to determine the parameter u in the problem. One simply demands that this collinear 
singularity vanishes: 



0, 



We are even forced to proceed this way to ensure the solubility of the problem. One 
has to fix UJ to the expression 



m. 



k 



j_ ml + k\ 

- H 



x(l — x) x'{l — x') 



(2.10) 



which indeed has the form of a kinetic energy. Although it can be shown that the 
two kinetic energies of this sum must be the same, I wrote T* here in the equivalent 
suggestive form of an average. It cannot be overstated that exactly this form of the 
expression ( p.lQ ) follows from the structure of the effective theory pO| . 



The physical interpretation of this procedure is not easy. One can imagine that the 
Q-space contains all interactions of the higher Fock sectors by effectively summing them. 
In this complicated effective interaction, we disregard all non-diagonal contributions. 
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But on the other hand, this is partly compensated by the special choice u = T* . One 
can therefore state that T* contains an approximation of the summed interactions of 
the higher Fock states. 



2.2 The effective integral equation 

We can proceed by calculating the matrix elements of the now well-defined effective 
Hamiltonian. It operates only in P-space. To read off the actual definition of the matrix 
elements, we have to write down the integral equation in which they are contained. The 
continuum version of Eq. (|2.1| ) is 



xil — x) 



- M„ ilJn{,x, k±; Ai, A2) 



r dx'd^k'^ {x, kr, Ai, X2\j{leYj{le)\x', fcj; A;, A^^) , f , D v^n 

/ . / 1 /,2 I i2\ I , WAX , A^, A2) — 



.ft^ Jd \ {II + ID ^xx'{l-x){l-x') 



(2.11) 

This form of the effective integral equation is very useful for a comparison with the ma- 
trix elements calculated in |^ , see also Appx. |^. However, a main topic of the present 
work is the investigation of rotational invariance, or in general Poincare invariance, of 
QED on the light-cone. In this context, it is helpful to write the interaction term of 
Eq. (|2.11| ) a covariant way. The integrand reads now 

9^ Sr f dx'(Pk'j_ j^(/e, K)jfi{le, ^e) j_ / / r, . / n , „n 

\'^,\'^''^ J XX'{1 — X){1 — X') 'e'e,/i 

which makes it obvious that the effective interaction^ 

j^(^e, Ae)j^(/e, Ae) 



eff 



is gauge invariant and a Lorentz scalar. We restrict the integration domain D using 
the covariant cutoff of Brodsky and Lepage @: 

^^^^<A^ + 4m^ (2.13) 
x[l — X) 

which allows for states having a kinetic energy below the cutoff A. The matrix elements 
are given in some detail in Appendix 0. The explicit functions are listed in Appendix 



■^Instead of providing the arguments of the currents by bras and kets hke in Eq. (2.11), we wrote 
them into brackets. 
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Figure 2.4: Diverging diagrams in P-space: self energy diagram. 





Figure 2.5: Diverging diagrams in P-space: contraction contributions. 



So far, we have considered only the problems that occur due to the projection of 
Q-space onto P-space. But already in P-space we have to face severe singularities stem- 
ming from the graphs of the electromagnetic self energy, Fig. ( |2.4| ), and the contraction 
graphs. Fig. 

A renormalization scheme is necessary to remove the divergencies. It is described 
in Appendix |^ and can be summarized as follows. Although each of these graphs is 
quadratically divergent in the cutoff A, their sum is only logarithmically divergent. The 
arguments used in the proof assume the formulation of the theory in the continuum. 
In the discretized theory the arguments do not necessarily hold. 

Before addressing to solve the effective integral equation (|2.11|) , let us briefly re- 
view the steps that lead to its derivation and compare the method of this work to 
other attempts to solve for the positronium spectrum. First, we restricted the Fock 
space to the sectors |ee) and leej). By this, we arrived at a positronium model with 
one dynamical photon. This constitutes a matrix equation (2x2 block Hamiltonian) . 
The advantage of this matrix equation is that all interactions in the Q-space are ex- 
plicitly taken care of, whereas in our model their contributions are occluded by the 
determination of the redundant parameter u. However, efforts to extract the spectrum 
and wavefunction from the matrix equation have failed. In the case of QED, Tang 
44 1 and Kaluza produced results with no clear significance. 



and Kaluza produced results with no clear significance. This was due to 
the convergence problems, because they did not include counterterms for the Coulomb 
singularity. The counterterm technique used in the present work is described in Ap- 
pendix Calculating the counterterms, one faces in general numerical integrations in 
two or three dimensions, which must be performed for each diagonal matrix element. 
The numerical effort increases tremendously with the matrix dimensions [E^ • We wish 
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to avoid the disadvantages of the matrix equation. Therefore we have performed the 
continuum hmit of the DLCQ formahsm and obtained a coupled system of integral 
equations. The effective integral equation ( p. 11 ) was derived by casting the effects of 
higher Fock states into the Q-space by fixing u to T*. Then the Q-space is projected 
onto the P-space as a second step. 

If one considers the non-relativistic limit of the interaction in the effective integral 
equation (|2.12| ), one arrives at the so-called light-cone Schrddinger equation. The main 
point in this calculation is the relation between the longitudinal momentum fraction x 
in light-cone coordinates and the equal-time momentum in z-direction, kz 



X 



x(^kz 




m? + k'^, + kl, 




The calculation of the non-relativistic limit is then straightforward. The solution of 
the light-cone Schrodinger equation, which has analytically integrable Coulomb coun- 



terterms, was given in ||2^, Fig. 8], shows very good convergence, and yields the correct 
eigenvalue spectrum. If one fully reduces the effective integral equation to the non- 
relativistic limit, one arrives at the well-known Coulomb equation in momentum space. 



solved numerically and discussed in Ref. p3 . 



2.3 The positronium mass spectrum 



The solution to Eq. (|2.11| ) was given in In contrast to the light-cone Schrodinger 



or the Coulomb equation, the counterterms for the Coulomb singularity cannot be 
calculated analytically. As can be seen in the helicity tables in Appendix 0, in the 
case Jz=^ essentially two different diagonal matrix elements occur: one for parallel. 



the other for anti-parallel helicities. Krautgartner |73] was able to integrate out 
analytically one of two variables contained in the continuous part of the counterterm 
for anti-parallel helicities. As a result, he only had to integrate over one dimension 
numerically. However, he did not succeed in analytically integrating out this variable 
from both matrix elements and consequently had to use the same counterterm for both 
diagonal matrix elements. The convergence and the spectra he obtained with this 
method were reasonably good. Indeed, the use of identical counterterms in this case 
(J2=0) is well justified since both functions have the same singularity structure and 
comparable values. We note that this becomes problematic in the case of non-vanishing 
Jx, where one has four distinct diagonal matrix elements, one of which is much smaller 
than the others. 

We calculate in this chapter the spectrum of the positronium model described in 
with an improved counterterm technique. This means a rigorous calculation of all 



four counterterms corresponding to the individual diagonal matrix elements. The prize 
is an entirely numerical two dimensional integration rather than a one dimensional 
analytic integration over the variable]^ cos 6', followed by a numerical integration of 

''For the definition of the variables used in the calculations on the computer, see Appx. |^. 
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the off-shell mass /i. The effect is an even better convergence of the spectra with the 
number of integration points, as one can see by comparing Fig. ( p.6|) with Fig. ( p.7|) . 
In particular, one notes a better convergence for a principal quantum number n>2 and 
a number of integration points A^<13. The lowest states for n=2 converge much better 
when the entirely numerically integrated counterterms are used. This is the result of the 
afore mentioned small, but distinct difference between the diagonal matrix elements. 
If only one counterterm, adjusted to the singularity structure of one special diagonal 
matrix element, is used, the Coulomb singularity of the other diagonal element is over- 
compensated. As a consequence, the calculated eigenvalue is smaller than it should be. 
This is exactly the effect observed in the spectra. 

We recall the analytic results^ for the singlet and triplet states @ll7T|, p. 10] to 
order O(a^) and write it in the form of Bethe and Salpeter ||7|, §23]: 



En,i = -\^y 



V? 32n4^V''^'-^ 21 ^\) n 



with the principal quantum number n and the Rydberg constant Ry = mj«^/2. The 
singlet terms have 

Q,5=o,J = 0, 



and the triplets 



31+4 



7 1-5,0 ('+i)p+3) 



if J = 1 + 1 



L l(2l~i) n — t -L- 

For a comparison of our results to "experiment", i.e. to perturbation theory with a 
strong coupling constant, a = 0.3, we have compiled the positronium mass spectrum 
in Table ( p. 21) . Usually one classifies the states according to their quantum numbers of 



^The "state of the art" theoretical results are given by Gupta et AL.[^. For the triplets they 
have 



J-Ry a3 



1 

An 

H — Ry a ^ ma 
6 n"* 



7 109 2, „ ^, _i 16, , , , 
a„ \ — ln2 + 61na InKo(n) 



15 



with the Salpeter terms 

9 

ais = -21n2-3, a2s = 

and the Bethe logarithms 

In fco(l) = 2.9841285558..., In fco(2) = 2.8117698931 ... , lnfco(3) 
A review of experimental results can be found in |70| Chap. 15]. 
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total angular momentum J, orbit angular momentum L, and total spin S. This is valid 
only for rotationally invariant systems, or, in our case, in the non-relativistic limit. 
These quantum numbers constitute the spectroscopic notation '^^'^^Lj. We choose a 
convention for the time reversal operation 

Vn\J,J.):={-lY'-'^^\J,J,) (2.14) 

to classify the states likewise. The singlet states are known as parapositronium, the 
triplet states as orthopositronium. We display the non-relativistic notation for the 
states in Table ( p.2| ) to make the comparison to other data easier. The eigenvalues are 
listed in the form of binding coefficients B^, defined as 

B„:^4^ (2.15) 

for all eigenvalues. 

The finite N error estimates given in the table were obtained by comparing the 
results for the maximum number of integration points with those for the next highest 
number of points. The actual errors are definitely smaller, because the eigenvalues con- 
verge exponentially with the number of integration points A^. We will comment on this 
in detail, when we have completed our model by introducing the annihilation channel 
in Chapter |. A word seems in order on the magnitude of the errors. Surprisingly, the 
largest errors are those of the states with the largest binding coefficients, in particular 
of the ground state. The explanation is that we work in momentum space. Conse- 
quently, the higher excited states, widely spread in coordinate space, are condensed in 
momentum space and therefore in the region of many integration points. 

Note the good agreement in Table (PT^), including excited states. The singlets 
^5*0 tend to be calculated as too weakly bound, especially for n=l and n=2. This 
effect is reversed for the triplets ^5*1. In principle, the whole bound state spectrum is 
accessible with our method, not only the first few states. To support this statement, it 
is instructive to investigate the properties of the higher states. 

The main features of a detailed study of the multiplet structure of the spectrum are 
the following and can be seen from Table 



There is a well-defined number of states for each fixed principal quantum number 
n. In the case considered here (7^=0), there are 4(n — 1) -f2 states. This number 
should be reproduced for as large an n as possible. It turns out that the multiplets 
contain the correct number of states in our model up to at least n = 5. 

Each state has defined quantum numbers concerning the charge conjugation ttc 
and T-parity ti-h, cf. Eqs. (|G.6D and ( p.7|) . These quantum numbers can be 



obtained from the non-relativistic notation "^^^^Lj by 

= (-1)^+^ and Tin = (-l)^+''+\ 



using the convention ( |2.14| ) . It is an important result that only those combinations 



of quantum numbers ttc and t^-h which are expected from the theory, occur in each 
set of eigenvalues for any given n in our calculations. 
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Figure 2.6: The spectrum of the effective integral equation for a = 0.3, A = 1.0 m/, = 0. The 

mass squared eigenvalues M,^ in units of the electron mass m'j are shown as functions of the number of 
integration points N = Ni = N2. The calculation was done using the entirely numerically integrated 
Coulomb counterterms. Note the improved convergence of the states with n>2 in comparison with 
Fig. (p.7\). The 100 lowest eigenvalues are displayed. 



The ordering of the muhiplets seems to have minor errors. For instance, the 2^5*0 
state and the 2^Pq state are permuted, cf. Table (|2.2|) : the S-state should be the 
lowest according to perturbation theory up to order 0{a^). This finite cutoff 
effect is explained in next paragraph. 



Cutoff dependence 

It is stated in Appendix ^ that there is, due to the renormalization scheme used, a 
logarithmic divergence of the eigenvalues with the cutoff A. This was already found 
in |]73|, Fig. 7]. In Fig. ( |2.8| ) and ( |2.9D the eigenvalues are shown as functions of the 
cutoff for n=l and n=2, respectively. It is obvious that the different eigenstates diverge 
linearly with the logarithm of the cutoff, and that the coefficients of these divergencies 
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Figure 2.7: The spectrum of the effective integral equation for a = 0.3, A = 1.0 m/, = 0. The 
mass squared eigenvalues Af^^ in units of the electron mass rrij are shown as functions of the number 
of integration points N = Ni ~ N2. The calculation was done using the half-analytical Coulomb 
counterterms of Ref. j2^. The convergence of states with n > 2 is not as good as in Fig. (2.6). The 
100 lowest eigenvalues are displayed. 



are different. It is reasonable to fit tlie curves of Fig. ( |2.8| ) witli a polynomial in log A, 
because this is the behavior expected from the renormalization scheme. If one omits 
the points for A > 20mj, because there the entirely numerical counterterm integrations 
become problematic, one gets good agreement with the calculated curves if one uses 

Msingieti^) = 3.90545 - 0.0350983 log A + 0.00745955 log^ A, 
Mtripieti^) = 3.90976 -0.01857871ogA + 0.00788614 log^ A. (2.16) 

The small coefficient of the log^ A term verifies the logarithmic dependence of the 
eigenvalues on the cutoff. We will see in Chapter | that the dependence on A becomes 
even weaker if one includes the annihilation channel. 

One notices several level crossings for n=2. As was stated in the last paragraph, 
the ordering of the eigenvalues of n=2 for A=1.0m/ turns out to be wrong for the 
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Cutoff: A 


Bs 


Bt 


Chf 


1 n 


1 04Q03Q64 


1 00046997 


1 34Q371 3 


X .o 


1 16373Q04 


1 06860Q34 


n 9fi494Q1 7 


o .u 


1 9^^70148 


1 1 m 1 1 S98 

X.XWXXXO.i<0 


n 49Q41 1 66 




1 9QQ7Rn^n 


1111 (S'^^TR 

X . X X X UOO 1 o 


n ^99(^9499 


7 9 


X . OZiC/'H: X C/ X 


1 1 1 7R97R9 

X . X X 1 OZi 1 OZi 




Q n 


1 '^^99'^QR9 


1 199'^'^(^^9 


u.uoouzuzo 


10.8 


1 37112216 


1 12596311 


68099735 


12.6 


1.38744792 


1.12904455 


0.71778713 


14.4 


1.40198469 


1.13175363 


0.75064183 


16.2 


1.41520247 


1.13419048 


0.78058886 


18.0 


1.42740143 


1.13641774 


0.80828803 


ETPT 


1.11812500 


0.99812500 


0.33333333 


0{a^ In a) 




0.23792985 



Table 2.1: The binding coefficients of the singlet (Bs) and the 
triplet states (Bt) for a = 0.3, A^i — 25, N2 — 21 as functions of 
the cutoff A in electron masses. Additionally the values for equal- 
time perturbation theory up to order 0{a'^) (ETPT) and up to 
order O(a^lna) (cf. Eq. |2.17|) are shown. 



two lowest states of this multiplet because of the crossing. That the levels do indeed 
cross can be proven by tracing them back to their sectors of definite C- and ?i-quantum 
numbers. A consequence of these crossings is the fact that the order of the eigenvalues 
is correct in the region between the crossings of the 2^S'o/2^Po and the 2^Pi/2^Po states, 
1.5 < A < 7. A further investigation of these crossings up to n=4 show, that the states 
with TTc= + 1 and 'n'n= — 1 are those that fall off fastest with A and tend to cross other 
levels. 

An important ratio of eigenvalues to be compared to results of other calculations is 
the hyperfine splitting. The hyperfine coefficient 



c 



hf 



triplet 

2 

3 + 



-M, 
1 
2 



singlet) / 

a 



TT 



ln2 + l^) 
9 ) 



12 



aMn a + Ka^ + K'a^ 



(2.17) 



was introduced by Fermi |7^ in 1930. Fermi calculated Chf = \ for hydrogen-like 
atomsQ, which is the exact result up to order O(a^). The second line of Eq. ( |2.17| ) 
shows the 'state of the art' result of equal time perturbation theory [|75|, p. 759]. The 
term in brackets is the contribution from the one-photon annihilation. This coefficient 
is listed in Table ( p.l|) together with the binding coefficients (|2.15|) . They will be com- 
pared to the values obtained by including the annihilation channel in Chapter ^. The 

^ He investigated the cases of Sodium and Caesium, e.g. he naturally did not consider the possible 
annihilation of a bound particle-antiparticle system. 
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Figure 2.8: Cutoff dependence of tlie triplet (up- 
per curve) and singlet (lower curve) ground state, 
a = 0.3. The cutoff is given in units of the elec- 
tron mass. 
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Figure 2.9: Cutoff dependence of the first ex- 
cited states (n=2) for a = 0.3. The cutoff is 
given in units of the electron mass. Note the level 
crossings. 



coefficients were calculated with the improved counterterms, as opposed to Table 
V.]. Comparing with the values given there, one notices that the singlet falls off slower 
with A for the old counterterms. Consequently the values for the hyperfine coefficient 
Chf are smaller there. The value of this coefficient, using Eq. ( p.l7|) without the annihi- 
lation contribution, is displayed, too. In the region of large couplings considered in this 
work, also higher orders in the coupling constant are important. Note the remarkable 
effect: the value Chf is 40% larger up to 0{a'^) than up to (^(a^lna)! 

Concluding, one can state that the best values as compared to equal time pertur- 
bation theory are obtained for A ~ 1.5 m^: the hyperfine splitting has the right order 
of magnitude, the order of the eigenvalues is correct, and the ground state is at the 
value of perturbative calculations. Although one can think of fitting the obtained data 
to the results of perturbation theory, I decided to follow the renormalization scheme of 
Appx. ^ to show the consistency of the positronium model described in this chapter, 
rather than to produce results competing with the accuracy of elaborate perturbative 
calculations. Moreover, it is not clear in the regime of a large coupling, how reliable 
the results of ETPT are, as we have seen from the significant dependence of the value 
for Chf on the order of perturbative calculations. 
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1 -L 


06268687 


06386407 ± 00030333 


—0 001177 


1.88 
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41 A 


-1 


-1 


0.06266009 


0.06314142 zb 0.00009618 


-0.000481 


0.77 


30 


43^2 


+ 1 


+ 1 


0.06266009 


0.06329675 zb 0.00011207 


-0.000637 


1.02 
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-1 


0.06267683 


0.06323357 ± 0.00011899 


-0.000557 


0.89 
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43^4 


+ 1 


+ 1 


0.06258754 


0.06309753 ± 0.00010310 


-0.000510 


0.81 
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5^5*0 


+ 1 


-1 


0.04134100 


0.04325281 ± 0.00125908 


-0.001912 


4.62 
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53^1 


-1 


+ 1 


0.04038100 


0.04291300 ± 0.00127772 


-0.002532 
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5IP1 


-1 


-1 


0.04038100 


0.04283966 ± 0.00064062 


-0.002459 


6.09 



Table 2.2: The positronium spectrum for a = 0.3, A = 1.0 m/, Ni = N2 = 21. The n'^^^'^Lj notation, 
the quantum numbers under charge conjugation, C, and T-parity, H, are shown in the first columns. 
The first row of binding coefficients (Betpt.u) comes from equal time perturbation theory calculations 
up to order 0{a'^). In the following row our results are listed with an estimate of finite N errors. Also 
shown is the difference AB := BETPT,n — Btheor,n- The last row contains the relative discrepancy to 
perturbation theory in percent. 
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Figure 2.10: The singlet wavefunction for anti-parallel spins as a function of the longitudinal mo- 
mentum fraction x and the transverse momentum k±, omitting the dependence on the angle Lp. The 
calculation was done with a = 0.3, A — l.OnifjJz —0,Ni= 41,iV2 = 11. 



2.4 Wavefunctions 

A big advantage of the Hamiltonian method apphed in the calculations of the spectrum 
is the fact that the wavefunctions of positronium are obtained in the same calculation 
as the spectrum. The wavefunctions for the singlet and the triplet components are 
displayed in this section. Similar plots in Ref. |73[ seem to indicate numerical problems 



of that work, because they show some internal structure. It turns out that this is due 
merely to mistakes in the graphing package used. The smoothing functions are quite 
sensitive to the data and boundary conditions employed^ 

The wavefunctions have two components: one with parallel (|t) and the other with 
anti-parallel (||) helicities. For the displayed singlet and triplet wavefunctions, the 
decomposition in terms of spin components (cf. Eqs. ||G.6|] and [ p.7|| ) read 



1 Ni {N+l)/2 

IV'sing(n)) = 2 E E [#I(T)4(i) -^*&I(i)4(T) - #ki)4(T) + ^*4(T)4(i)] |o), 



=1 j=i 

Ni {N+l)/2 



-I iVi yiy^L//^ 

IV'sing(TT)) = E [#I(T)4(T) + ^*fel(i)4(i) -#kT)4(T) -^*4(i)4(i)] |0) 



2 i=i j=i 



For our purpose GLE 3.3 by C. Pugmire yielded the best results together with the smoothing 
function grid of PLOTDATA. 
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Figure 2.11: The singlet wavefunction for parallel spins as a function of the longitudinal momentum 
fraction x and the transverse momentum fc^, omitting the dependence on the angle Lp. The calculation 
was done with a = 0.3, A = 1.0 m/, 0, iVi Al,N2 = 11. 



. ATi {N+l)/2 

Iv^tnp(Ti)) = ^ E E [#i(T)4(i) + ^*&I(i)4(T) + #ki)4(T) + ^*4(T)4(i)] |o), 



-| ivi v-i'^J-;/^ 

IV'tnp(TT)) = E E [#I(T)4(T) - rb\{[)d\{i) + #kT)4(T) - rhl{[)d\{[)\ |0), 

The wavefunctions are normahzed in the polar coordinates /x and cos 9, such that 



Nx (iV2+l)/2 

E E [^(/ij,cos6'j; Ai = A2) + V'(/ii,cos6'j; Ai = -A2)] = 1. (2.18) 
i=i j=i 



where 



61 = /c_L, = cos6', yj), 

(i2 = (i^ (1 — x, /c_L, = — cos 6*, + tt). 

and discretized variables /ij, cos^^j are used. For all details on the numeric aspects of 
the calculations see Appx. ^ 

Figure ( p^.l4D shows the decrease of the singlet wavefunction with anti-parallel he- 
licities as the off-shell mass increases. The graph, calculated with entirely numerical 
counterterms, is almost the same as that in Fig. (11), bottom] with half-analytical 
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Figure 2.12: The triplet wavefunction for anti-parallel spins as a function of the longitudinal mo- 
mentum fraction x and the transverse momentum fc^, omitting the dependence on the angle Lp. The 
calculation was done with a = 0.3, A ~ 1.0mf,Jz ~Q,Ni^ 41,iV2 ~ 11. 




Figure 2.13: The triplet wavefunction for parallel spins as a function of the longitudinal momentum 
fraction x and the transverse momentum kj_, omitting the dependence on the angle if. The calculation 
was done with a = 0.3, A = 1.0 m/, = 0, TVi = 41,iV2 = 11. 
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Figure 2.14: The decrease of the Jz — singlet ground state with anti-parallel helicities as function of 
the off-shell mass variable fi, Eq. mu. The parameters used are Ni = 41, N2 = 11, A = 1.0 m/, a = 0.3. 
Plotted are the six different decreases corresponding to the six non-positive values of the discretized 
angle variable cos6. One notices the deviations from rotational symmetry for /i > 10. 



counterterms, since a large number of integration points (A^^i = 41, N2 = 11) were used 
to calculate the results. It can be seen that rotational invariance is broken. This com- 
ponent of the S-wave state does not depend on any angle and should therefore decrease 
independently of cos^. The broken rotational invariance is no surprise. One reason 
is the fact that the associated operators Fi and F2 of rotations around the transverse 
axes are dynamical, i.e. contain the interaction, cf. Table ( |3.1| ). The other argument 
uses the transformation to equal-time coordinates 

^^(i-I)a-. (2.19) 

One sees directly that the term proportional to -p- breaks rotational invariance. 

It is worth mentioning that the breaking of rotational invariance is noticeable only 
for large cutoffs A. With a cutoff of A = 1.0 m^, the deviation of curves with different 
cos^^ is not visible. In the next chapter, we investigate the same properties of the 
corresponding wavefunction for 7^=1 . 



dxd k I 



dx 
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Chapter 3 

Angular momentum in front form 
dynamics 



3.1 The Poincare group on the hght-cone 

To provide the theoretical background for our numerical calculations, we review the 
properties of angular momenta on the light-cone. The material presented here was 



inspired by a section on the same subject in |76 



3.1.1 Unitary representations of the Poincare group 

The Poincare group has 10 generators: the four-momentum P'^, and the generalized 
angular momenta, constituting an anti-symmetric tensor M'^'^ = —M'^^. They are 
subject to the commutation relations 

[P^,P''] = 0, 

The Poincare group has two Casimir operators, i.e. two invariant functions of the 
generators. One is the mass (squared) operator 



:= P^P, 



A" 



which cannot be confused with the mc^exec? tensor of angular momenta M^^ . The other 
is the square of the Pauli-Lubanski vector 

1 

— ( 
2 



:= -e^^^^PPM^'r 



The components of this vector obey the commutation relations 

[WP,W']=ie,,,„WPP% (3.1) 
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and commute with P^^. The spin operator j can be defined ||76| , Eq. (6)] as 

j-:=<(P)^, ^ = 1,2,3. (3.2) 

It is therefore a hnear function of the Pauli-Lubanski vector. The coefficients Ui are 
three orthonormal, operator valued, space-like basis vectors, orthogonal to P^. They 
have the properties 

<(P)«,,,(P) = and f:<(P)<(P) = g^-" - (3.3) 

i=i 

The basis vectors transform under a Lorentz transformation like 

f/t(A)M,(P)f/(A) = Ui{AP) ^ Au,{P). 

These two basis systems are related by a rotation, since Mj(AP) and Aui{P) are both 
orthonormal space-like vectors, orthogonal to AP'^. This transformation is referred to 
as a Wigner rotation by Coester W^ 



One can show that the Lorentz transformation of the spin operator is equivalent to a 
Wigner rotation of its components 

U\A)jiP)UiA)=nw{A,P)3. 



This follows from the definition ( |3.2| ) of the spin operator and holds for every vector 
6j(P) = -uf (P)V^, where is an arbitrary four- vector. The spin operator has the 
correct SU(2) commutation relations, as can be proven using ( |3.1| ) and the properties 
of the basis vectors, Eqs. (|3.3| ). Lorentz transformations A have the general property 



A/ g,^A,- = g,,. (3.4) 
If one defines Uq := one can re- write Eq. (|3.3| ) into a four dimensional equation 



which has the same structure as Eq. ( |3.4|) . The operator valued matrix is thus an 
50(1, 3) representation of some Lorentz transformation B{P) which can be interpreted 
with some care as the Lorentz transformation to the rest frame of P: u^{P)Pi, = g^^M. 
The spin operator Jean be written, apart from the form given in Table ( |3.1|) , in terms 
of the basis vectors 

We are interested in bound state calculations, i.e. in composite systems. One can define 
the relative momenta kn of a subsystem by 

/c„, := <(P)p„,„ p'^^ = Y.u>t{P)K,i, 

i 

where P = J2nPn is the (absolute) momentum of the system. 
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3.1.2 Connection between canonical spin and front form he- 
licity 

So far, everything was independent of the choice of the basis vectors Ua{P), i-e. of the 
form of the Hamihonian dynamics. The imphcations for instant form and front form, 
the two kinds of dynamics of interest for the present work, are compared now. The 
major features of the two types of dynamics are hsted in Table (|3.1|). 

The instant form, i.e. canonical, choice of the basis Ua{P) is related to the front 
form basis Ua{P) by the so-called Melosh rotation [7S] 

nMiP)^r■=U^iP)■UJiP). 

Its SU{2) representation is 

M + P+ -ia -{zx Pi_) 



D^/^[nM{P)] 



(M + P+)2 + Pi 



The essential observation for a relation between instant form spin and light-cone helicity 
is that the canonical spin of the n-th particle of a composite system can be expressed 
in the front form basis 

Sn,^ = U^iK) ■ B{P)^ = TlMiPMh) " B{P)^. 

IVln -'"n 

The translation of the total angular momentum from equal time to light-cone coordi- 
nates is therefore obtained by a Melosh rotation, cf. Table (p . 1|) . 



3.1.3 Field theories on the light-cone 

There are two different ways to construct front form particle dynamics. One approach 
starts from the mass and spin operators which have to fulfill some constraints on their 
commutation relations. The three Hamiltonians of front form dynamics can then be 
expressed as functions of the kinematic generators and of the mass and spin operators 



M^ + Pl 
P+ 

_ 2 (Mj, + P,J3) P- 2P, 
^1 - + J^Ei + —K,, 



P 



_ 2 {Mh + Pih) . P p ,2P2 

^2 + + 



The dynamical operators commute with each other 

[P-,FJ = [F1,F2]=0. (3.5) 



31 



Chapter 3. Angular momentum in front form dynamics 



The other approach is that of Fock space field theories. Here, the fundamental quan- 
tities are the Hamiltonians which are derived from a Lagrangian density. The spin 
operators are then formulated in terms of the Hamiltonians 

, - M - ^■(^±y< P±) 

Jz ^VIz , 

■'^ ~ ~W ~ MP+ ' 

where 

W>± = ^P^z X (f - + i (P+ - P-) zxE^-[zx Pj_) K,, 

and z is the unit vector in z-direction. Very important is the fact that Poincare invari- 
ance is destroyed, as soon as truncations of the Fock space or regularizations of Fock 
sectors are implemented. In particular, the requirement of an invariance under rota- 
tions around the transverse axes is difficult to fulfill since the corresponding operators 
are complicated and involve the interaction. Even worse, in a truncated Fock space 
formalism the full Poincare invariance is absent if it is not restored by an additional 
[ad hoc) prescription ||76|, p. 11]! 

It is therefore impossible for our positronium model to be Poincare invariant. I shall 
rather show its covariance by looking at the results, i. e. at the physical observables such 
as the invariant mass spectrum. If full rotational invariance is restored in the solution, 
the states of same total angular momentum J but different Jz, become degenerate. 

The most direct way would of course be to construct the operators Fi and F2 ex- 
plicitly and the to diagonalize the operator of total angular momentum J'^. This has 
not been done up to now. Because of the vanishing commutators of the dynamical op- 
erators with P~, Eq. ( p.5| ), it is clear that the diagonalization of the rotation operators 
will be much simpler with an already diagonal Hamiltonian Hlc'- only the states with 
degenerate mass eigenvalues will be coupled by F±. 

However, we restrict ourselves in the present work to the calculation of the spectrum 
of the light-cone Hamiltonian and will ffnd that we can classify the eigenstates with 
regard to even without constructing or diagonalizing the rotation operators F^. 
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instant form dynamics 


front form dynamics 


Kinematic 
operators 


6 operators: 

P, J := {M^^,M^\M^^) 


7 operators: 

P,K,:^Im+-,E^:=M^+,J, 


Dynamic 
operators 


4 operators: 

P^,K 


3 operators: 

p-,Px := M^- 


Conditions 
on basis 


B[P) are rotationless Lorentz 
transformations 


B{P) are null-plane boosts 


Basis 




u± :^ e±+ n, 
ni^Pti 

M ( (n^'P.)P\ 


Wigner 
rotation 


nw{n,p) =n 




Angular 
momentum 


^ W PWo 
~ M M{M + po) 


j^n^{p)j 


Spin 
operator 


Si = 1Zw{P:Pi)ji 




Notations 


e/:=g/, n^:={l,-z), Aboost ^ {E^, K,} 



Table 3.1: Major properties of instant form and front form dynamics. 
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3.2 The Hamiltonian matrix with general 



The material presented in the previous section shows that the definition of angular 
momentum operators in front form dynamics is problematic since they include the in- 
teraction. It is therefore a subject of its own merit to study the properties of angular 
momentum within a well-defined model on the light-cone. Consequently, we investi- 
gate the case of a non-vanishing 2;-component of the total angular momentum in our 
positronium model. 

The way to proceed is inferred from the definition of the integral used to integrate 
out the angular degree of freedom (ip) and substitute it with the discrete quantum 
number J., 



{x,k±; J^, Ai, A2|Kfr|a;', A;^; J^,A'^,A2):= 



277 Jo 



2tt 



(3.6) 



dip e 



/ dip' e'^^"^ {x,k±,(p; Ai, A2I Ksla;', ^1, v?'; A'i,A2). 
Jo 



We listed the general matrix elements {x,ip] Si, S2\Ves\x' ,ip'; s[, S2) in Table ( [l*'.!! ). In 
general, the functions displayed there contain a dependence on the angles. Hence, it 
is not clear how the general ip dependence of the matrix elements will look like, if one 
inserts an arbitrary = Jz — into Eq. (|3.6|) . Fortunately, a simple scheme can be 
set up to construct the functions for all = n, n G Z. In particular, one can prove that 
the matrix elements can only depend on the difference ip — ip' . The general function 
has the shape 



\ p2iv /■27r 

/c_L,x', fc^; Ai, A2) = — / dip j dip 
An Jo Jo 



^ J F{ x, k±,x',k'_^]Xi,\2) J, 
- 2k±k'j_ cos(v9 — ip') 



where n G Z and 



(X — X 



(X — X 

2^ 



1 

\xx' 

k'? 



+ 



;i-x)(l-a;') 



+ ki +k\ 



-k' 



1 



1 



X 



.1 — x' x' 

It is straightforward to evaluate this expression with the decomposition exp{±zx} 
cos X ±i sin x and the integrals 



2lT 



2-E 



27r Jo 

1 

2^ 

Here, the definitions 



d^' 



2tv 



dip / dip' 



A 



cos{n{ip — ip')} 
a — 2k±k'j_ cos((/9 — ip') 

sm{n{ip — ip')} 
a — 2k±k'j_ cos(v9 — ip') 

and 



27r(-A)-l"l+i 
0. 



B 



kik', 



n\ 



4A;2 k'? 



B 



1 



1 - aA) 



were used. Using these relations, one can calculate the matrix elements for arbitrary 
X. The results are listed in Table ('IF.2D. 
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3.3 The positronium spectrum for general 

The positronium spectrum is calculated numerically from the matrix elements of Table 
([F.2|) . Although we implemented a more general counterterm technology, described in 
Chapter 1, which automatically accounts for the new features of the diagonal matrix 
elements for 7^ 0, we must be even more careful. In the former calculations (J^ = 0), 
we made use of the discrete symmetries C and as described in Appendix §2. These 
symmetry properties are not explicitly conserved for the more general case 7^ 0. 
To be careful, we ignore possible point symmetries in the problem, and solve for the 
spectrum without symmetrizing the Hamiltonian. We will see in the next section how 
justified our extra care was. The numerical effort increases enormously. With the 
unsymmetrized Hamiltonian, the dimensions of the matrices to be diagonalized are 
four times larger. Since the number of operations grows with the third power of the 
matrix dimensions in a standard diagonalization algorithm, the CPU time used is 16 
times longer. 

We calculated spectra for the seven different values, = — 3, — 2, . . . , +3. It is 
found that the eigenvalues are identical for J^= + n and Jz= — n a.s can be seen from 
Fig. The individual spectra and the convergence can be seen from Figures (|3.2|) - 

(|3.4|) . In particular, one notices that the singlet ground state is absent from all three 
plots, and that for 7^=2 and 7^=3 even more states are absent. The explanation is 
the impossibility to have states with = n in multiplets with J < n. The numerical 
stability {i.e. convergence) is very good: in each of the figures the lowest eigenvalue is 
almost independent of the number of integration points. In fact, the eigenvalues con- 
verge exponentially with the number of integration points, as we will show in Chapter 
|. This is the more surprising, as we adjusted the Coulomb countertermQ, based on 
the non-relativistic ground state wavefunction, and the excited state wavefunctions are 
quite different. 

Fig. ( |3.1| ), the summary of the spectra for different J^, is central to this chapter. 
What can we learn from it? It has two prominent features. Firstly, there are multi- 
plets of states with different Jz that are degenerate. We shall discuss the numerical 
evidence in what follows. Secondly, there is a limited odd number of degenerate states 
for each eigenvalue. The interpretation of these facts is obvious. The positronium mass 
spectrum is a physical observable, Lorentz invariant and therefore independent of the 
mathematical algorithm applied and of the Lorentz frame used. Central forces are ro- 
tationally invariant, and this should be observed in the spectrum of an electromagnetic 
bound system, too. Rotational symmetry tells us that there has to be a defined number 
of degenerate states for each fixed value of the total angular momentum J. Conversely, 
since this is exactly what we observe here, we can infer the quantum number J from 
the number of degenerate states for a fixed eigenvalue M^. Concluding, the 1, 3, 5, . . . 
degenerate states constitute the singlets, triplets, pentuplets, . . . of a J = 0, 1, 2, . . . 
multiplet. 

In Table ( p.2|) , the spectrum obtained for J^=l is compiled to compare it with the 
^i.e., the converging function g{p,p') was chosen in this manner, cf. Appendix 
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Figure 3.1: Compiled spectra of positronium with different Jz = —3, —2, . . . , +3. All spectra have 
been calculated with a = 0.3, A = 1.0m/,iVi = N2 = 21. The mass squared eigenvalues in units 

of the electron mass mj- arc shown. The notation for the states is The resolution of the plot 

is inadequate for the multiplcts. Nevertheless, the numerical degeneracy of the three triplet ground 
states '^Sj^^,'^Si, and ■^Sl becomes very clear. 
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Figure 3.2: The spectrum of the effective integral equation for a — 0.3, A = 1.0 m/, Jz = 1. The 
mass squared eigenvalues M,j in units of the electron mass are shown as functions of the number of 
integration points N = Ni = N2- The calculation was done using the entirely numerically integrated 
Coulomb counterterms. Note that the singlet groun d st ate is absent. For n=l only the triplet '^S'l 
survives the projection on the Jz=l sector. Cf. Fig. ( |2.6| ). 




3.95 




3.95 



Figure 3.3: The spectrum of the effective inte- 
gral equation for a = 0.3, A = 1.0m/, = 2. 
Note the total absence of any n = 1 state. 



Figure 3.4: The spectrum of the effective in- 
tegral equation for a — 0.3, A = 1.0m/, Jz=3. 
Even more states are not compatible with Jz=3. 
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Figure 3.5: Deviation of corresponding eigenvalues for Jz=0 and Jz=l mul- 
tiplets {a = 0.3, A = 1.0 to/) with growing number of integration points. The 
graphs show AAP M^{Jz=0)- M^{Jz = l) for the states l^^i (A), 2^Pi (□), 
and 2^ Pi (o). 



eigenvalues for Jz=0. Apart from the absence of the states with J=0, the table displays 
an almost complete coincidence of the corresponding states. Only for the triplet 1^5*1, 
the gap between the two states of different is bigger, though the numerical error is 
actually smaller. 

We have to investigate the significance of the degeneracy with respect to the number 
of integration points and with respect to the cutoff A. To find out if this discrepancy 
is merely a numerical artifact, or a property of the positronium model, consider Figure 
(|3.5|). Here, the mass (squared) discrepancy between the Jz=0 and Jz=l eigenvalues 
is plotted versus the number of integration points N for three different states. One 
notices the convergence of AM^(l'^S'i) with N. The curve does not converge to zero, 
as one would want, but to a value of AM^(^S'i) ~ — 5 x 10^^. The mass gap for 
2^Pi does, however, go to zero as N grows large. The mass gap of the other helicity- 
triplet state 2^ Pi has the same increase as AM^(l^S'i), if one disregards the behavior 
of the corresponding graph for the untrustworthy values N=5 and N=7. It converges 
to 4 X 10"^ as N increases. We mention that Kaluza and Pirner find that in 
light-cone perturbation theory there is a discrepancy between the case of Jz=0 and 
7^=1. This is due to the perturbative method applied there and even expected from 
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Figure 3.6: Deviation of corresponding eigenvalues for Jj=0 and Jz=^ {a = 
0.3) as functions of the cutoff A. The graphs show AM^ :== M^iJ^^O) - 
Ml{Jz=l) for the states l^Si (A) [upper box], and 23Pi (□), and 2^Pi{o) 
[lower box]. Note that the scales of the two boxes differ by a factor of 100! 



the point of view taken in the present work. Perturbation theory to any finite order 
breaks the symmetries of the theory. Since rotations contain the interaction in front 
form dynamics, the associated symmetry will be broken in a perturbative approach. 

Cutoff dependence 

For the cutoff dependence of the eigenvalues for non-vanishing the statements given 
in the context of Chapter ^ hold equally. A main result of the present work is the 
documentation of the restoration of rotational symmetry in the solution in front form 
dynamics. Consequently, the investigation of the degeneracy of eigenvalues of same 
total angular momentum, but different Jz is crucial. The discrepancy between cor- 
responding eigenvalues of Jz=0 and J^=l as functions of the cutoff is displayed in 



Fig. (|3.6|) . One notices that there is a slight deviation of the cutoff dependence of the 
triplets l^Si for different Jz- The deviations of the other states (2^Pi and 2^Pi) are 
suppressed by roughly a factor of 100. 

This weak dependence of AM^ on the cutoff will be suppressed even more, if the 
annihilation channel is included. 
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32 


A^Fi 


0.063098 ± 0.000103 


0.063422 ± 0.000158 
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5^5*0 
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34 


53^1 
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0.042915 ± 0.000740 


-0.19 


35 


5^Pi 


0.042840 ± 0.000641 
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Table 3.2: The positronium spectrum for a = 0.3, A = 1.0mf,Ni = N2 ^ 21. The non-rclativistic 
notation for the terms and the binding coefScients for Jz=0 and Jz= + 1 are shown. The numerical 
errors are estimated from the difference between the values for maximum and next to maximum 
number of integration points. In fact they are smaller, because the eigenvalues converge exponentially 
with the number of integration points. 
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3.4 Wavefunctions 

The 7^ wavefunctions show more structure than those with Jz=0, due to the 
lower symmetry. In the Jz=0 case, there are essentially two different components of 
the wavefunctions: one for parallel, the other for anti-parallel helicities. Consequently, 
there are only two plots shown for the singlet and triplet wavefunctions in Chapter 1: 
Figs. ( CT ), (2.11) and Figs. (2.12), ( pTTSD , respectively. 

This is not a consequence of the symmetrized Hamiltonian. If the non-symmetrized 
Jz=0 Hamiltonian matrix is diagonalized, the same eigenf unctions are found as in the 
symmetric case, but twice as many occur. The only difference is a sign, depending on 
the parity quantum number. 

When Jz^O, we encounter four distinct components of each of the eigenfunctions, 
corresponding to four different helicity combinations. We elaborate on this subject by 
considering the components of the triplet wavefunction for Jz=^, Figs. (|3.8| )(a)-(d), and 
that of the next higher state. Figs. ( |3.9| )(a)-(d). In both cases, the components for anti- 
parallel helicities are identical, though displayed differently to show their full shape. For 
the triplet, the components with parallel helicities have nothing in common: the (tt)- 
component peaks at x=0.5, kj_=0 and is almost rotationally invariant, whereas the (ID- 
component vanishes at kj_=0 and is shaped more like the components with anti-parallel 
helicities. Note the extremely differing peak values: the anti-parallel components are 
suppressed by a factor of 40, compared to the (||)-component, the (||)-component 
even by a factor of 1400! As with the Jz=0 sector, Eq. ( |2.18| ), the normalization is 



Ni N2 

V^(/Xi,cos%; Ai,A2) = 1. 

i=l j=l Ai,A2 

The missing symmetry of the components is a consequence of the missing symmetry 
of the Jz sector with respect to T-parity. This property is found in all wavefunc- 
tions of these sectors. We have displayed here the wavefunction(s) of the next excited 
state to show another important fact: the wavefunction Fig. ( |3.9|) (a) has a small, but 
noticeable deviation from the reflection symmetry with respect to the x = 0.5 plane. 
This symmetry around x=0.5 is due to the charge conjugation symmetry C: if one 
permutes particle and antiparticle, one substitutes x with 1 — x. The fact that this 
property is respected by all wavefunctions other than Fig. ( |3.9| )(a) shows, that this 
symmetry is not broken, even not in the Jzy^O sector. The slight deviations can be 
explained as numerical errors, or more likely, as errors due to the grid function of the 
graphing package used. Some examples for higher excitations and larger are given 
in Fig.(p:TO|). 

The decrease of the triplet wavefunction l'^S'i(lt) with parallel helicities is plotted in 
Fig. ( |3.7| ) . These curves for different cos 9j have to be compared to those of Fig. ( |2.14| ) . 
In both cases, rotational symmetry is broken, since the decrease of the wavefunction 
is not isotropic but depends on the angle 6. There are some differences between the 
sectors Jz=0 and Jz=^- One is the fact that the smallest value of the wavefunction 
for 7^=1 is roughly 9 x 10~^, whereas for Jz=0 its approximately three times smaller. 
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Figure 3.7: The decrease of the Jz = +1 triplet ground state wavefunction 
with parallel helicities as a function of the momentum variable /i. The parame- 
ters are a = 0.3, A = 20. Om/, Ni = 41, N2 = 11. There are six different curves 
corresponding to six values of the discretized angle variable, 9: they show the 
decrease in V' with increasing jjL. Notice the deviations from rotational symme- 
try for II > 10. 



Another difference is the value of ix from that on the deviations in the decreases become 
noticeable. For Jz=^ this value is at /i ~ 10, contrary to ~ 3 in the Jz=l sector. 
Moreover, the curves of different cos 6* seem to be grouped for Jz=l. In any case, 
the important result is the same as in the case of J^^O: the wavefunctions are not 
rotationally invariant. 
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Figure 3.8: The triplet ground state wavefunction for Jj= + 1 as a function of the longitudinal 
momentum fraction x and the transverse momentum fc^, omitting the dependence on the angle (p. 
The calculation was done with a = 0.3, A = 1.0to/,7Vi = 41,iV2 = 11. Shown are: (a) (tT)-component, 
(b) (ti)-component, (c) (J,t)-component, (d) (J,J.)-component. 
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Figure 3.9: The 2^Pi wavcfunction for Jz=+1 as a function of tlie longitudinal momentum fraction x 
and the transverse momentum kj_, omitting the dependence on the angle ip. The calculation was done 
with a = 0.3, A = l.Orn/, A^i = A1,N2 = 11. Shown are: (a) (tt)-component, (b) ("f |,)-component, (c) 
(iT)-component, (d) (H)-component. 
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Figure 3.10: Components of wavefunctions for larger values of Jz as a function of the longitudinal 
momentum fraction x and the transverse momentum k±, omitting the dependence on the angle Lp. The 
calculation was done with a = 0.3, A = 1.0m/,iVi = 41,7V2 = H- Shown are: (a) ("f |)-component of 
2^Si , (b) (ii)-component of 235*1, (c) (it)-component of S^Dg, (d) (ii)-component of 3i£»2- 
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Chapter 4 

The annihilation channel 



4.1 Introduction 

So far, we have considered an effective Fock space, consisting of two sectors |ee) and 
1 667). We shall justify the absence of any higheiQ Fock state in Chapter ^ from the 



structure of the applied formalism of effective interactions. The general formalism 
is set up for a non-abelian SU (N) gauge theory. Unlike QED, the one boson state is 
absent in these theories because of color neutrality. Nevertheless, one has to take care 
of the one photon state I7) in QED, and the proceeding to include it is as follows. 

Firstly, it is important to notice not only the differences between the QED Table 
( |4.1| ) and the analogous QCD table in ||2y. Fig. 1], but also the similarities. Some of 
the graphs occuring in the QCD case are absent in QED, in particular those graphs 
with a three- or four-boson interaction and the instantaneous interactions connecting 
four bosons. But, although an additional sector occurs as a first row and a first column 
in the QED table, neither is a change in the higher Fock sectors observed, nor is the 
ordering altered in any way. 

In addition to Eqs. ( p.2|) and (|2.3|) , the A^-space {i.e. the sector containing the I7) 
states) is added to the system. The corresponding projector is 

N:= Y: l(7)n)((7)n|. 

n 

all QN 

The whole procedure of subsequent projections of higher Fock states onto the remaining 
Hilbert space of states can be carried out like before until one arrives at a (2 x 2) matrix 
analogous to (|2.4| ), but operating in the A^- and P-space rather than in the P- and Q- 
space 

From Table ([4.1| ) one can read off the interaction of the one photon state with all other 
sectors: the vertex interaction annihilates the photon into an electron-positron pair 

-'^Or rather the substitution of effects of higher Fock sectors with the use of effe ctive matrix elements 



in the remaining sectors. "Higher" here in the sense of ascending n in Table (|4.l| ) 
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Figure 4.1: The Hamiltonian matrix for QED. The sectors n are numbered starting at zero. The 
vertex, seagull and fork interactions are denoted by V, S, F respectively. Diagonal matrix elements 
are symbolized by •. This table is courtesy of H.-C. Pauli. 



and a fork interaction scatters it into the sector |ee7). The latter interaction is aheady 
contained in the effective interaction Eq. ([4.1| ) because of the projection of the Q-space. 

Although we always projected the higher Fock sectors onto the lower ones up to 
now, one is, of course, free to project the (lower) |7)-sector onto the (higher) |ee)-sector, 
and obtains 

Hcsi.^) = Hpp + HpN — — Hnp + HpQ — — Hqp 

Of course we do this for convenience; we could just as well solve the eigenvalue 
problem of Eq. ( |4.1| ). The projection is depicted in two graphs. One is the dynamic 
annihilation graph, Fig. (|4.2| ), the other is the corresponding seagull annihilation graph, 




1— X,— k^ 1— x',— k^' 1— X,— kx 1— x',— k' 



Figure 4.2: The dynamical annihilation Figure 4.3: The seagull graph of the annihi- 

graph. lation channel. 
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Fig. ( |4.3| ). The latter is a P-space graph and was absent before because of the gauge 
principle of Tang. 

4.2 Calculation of the matrix elements 

The calculation is straightforward. The matrix elements of the canonical Hamiltonian 
are given in Appendix 0. Up to now, we had to consider merely the first type of vertex 
interaction , cf. Table (|C.1|) , where a photon is irradiated from a fermion. Now the other 
graph must be evaluated. The calculation of the graphs [Fig. (|4.2| )-(|0|)] is performed 
in the (J^=0) and (J^= ± 1) sectors. The graph is absent in all other sectors because 
of the helicity of the photon: no angular momentum larger than J = 1 is possible. The 
functions derived depend on the light-cone momenta {x, k±). They are given in Table 
(|F.3|) and have to be added to those of Table ([F.2|) . The energy denominator in the 
one photon sector is simple because the photon has zero mass: 



where we used the definition ( p. 10 ). Note that this denominator does not depend on 



the directions of the vectors k'j_, i.e. on the angles if'. 

The calculation of the dynamic annihilation graph follows the steps described in 
Appendix 0. We mention only some aspects here. The transversal photon momentum 
vanishes: 

k±^ = /c_i_e + k±g = 0, 

and the longitudinal momentum is unity. Consequently, all expressions of the form 

vanish and Xj ^Jx\ = 1 becomes trivial. The matrix elements of the vertex interaction 
Vg^qq [Table ( |C2| )], split up into their three components with different helicity factors, 
are 

~ ^ ^ +Aii:+Ai 



{ee\VA\l) = + 

X Jb 

Jb 



The complete matrix elements are obtained by symmetry, since (ee|Vi|7) = (71^4 166)" 
To substitute the angles, we multiply the functions (ee|V^j|7) G{uj) (7|Vj|ee) by 
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final : initial 


(A'i,A'2) =TT 


(A'i,A'2) =n 


(A'i,A'2) =iT 


(Ai,A'2) =ii 


(Ai,A2) =TT 


Waa 


Wab 


Wac 





(Ai,A2) =Ti 


Wba 


Wbb 


Wbc 





(Ai,A2) =iT 


WcA 


WcB 


Wcc 





(Ai,A2) =ii 















Table 4.1: Symbolic helicity table for the dynami c an nihilation graph. The functions Wu 
are identical with the expressions Fi listed in Table ( F.3 ). Here, terms proportional to <5|j^|,o 
are omitted. 



following Eq. ( |3.6|) and integrate over ip and It turns out that because of the simple 
energy denominator (O), the only dependence on the angles comes from the factor 
(|4.3|) and is proportional to cos{(/? — or sinjyj — ^p'}. As a result, all matrix elements 
of the dynamic annihilation graph for 7^=0 vanish when integrated over the angles. 
Only for 7^= ± 1 do some of the matrix elements survive the integration. 

Because of the combination of matrix elements with the factor ( ^.3| ), two types of 
functions emerge for 3^= ± 1: one is independent of the angles, the other has a depen- 
dence proportional to exp{±2i((y9 — p')}. The latter vanishes after angular integration. 
The helicity table (^4.1|) is given to illustrate the helicity dependencies. It holds for 
Jz= + 1. The analogous table for Jz= — 1 is obtained by the operation 



Wi,{J,=+l- Ai,A2) = -\iW,,{J,= -l- -Ai,-A2). 

The simple kinematics {x^ + Xg = 1) of the seagull annihilation graph. Fig. 
in a constant contribution of this graph to the Hamiltonian matrix. It is 



(ee|S'|ee) - 

Because of its helicity factors, the graph acts only between states with 

S, = S'= 0. 



), result 
(4.4) 

(4.5) 



This means that the seagull graph does not contribute when ^ because it has a 
factor proportional to {ip — ip') resulting from ( [4.5| ). A rather surprising consequence is 
that the dynamic graph is the only annihilation channel contribution in the = ±1 
sector, whereas in the Jz = sector the instantaneous graph is the only one. This is 
necessary for rotational invariance: both diagrams must yield the same value, though 
one shows much more structure than the other, since degeneracy of the orthopositro- 
nium ground state with respect to Jz is found without the annihilation channel and 
inclusion must not destroy it. 

At first glance, there seems to be a manifest breaking of rotational symmetry: the 
helicity table ( |4.1| ) separates between states with (Ai,A2) = (tt) and (Ai,A2) = (ii). 
But this is only a consequence of the integration over the angles: for = +1 the (ID- 
combination gives no contribution, and likewise does the (||)-combination for Jz = —1. 
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Figure 4.4: The positronium spectrum including the annihilation channel: (a) Jz—0 sector, (b) 
Jz= + 1 sector. Parameters of the calculation: a — 0.3, A — 1.0 m/. The mass squared eigenvalues 
Afj^ in units of the electron mass are shown as functions of the number of integration points 
N = Ni = N2. The triplet states, especially l^Si, are lifted up, the singlet mass eigenvalues are the 
same as without annihilation channel. Cf. (E^) and Table (4.3) . 



4.3 Spectrum including the annihilation channel 

The spectrum including the annihilation channel shows the expected properties: the 
singlet eigenvalues remain unchanged, only the triplet states do change at all. An 
essential point in the actual calculations is that one has to use the same counterterms 
for the Coulomb singularity as used without the annihilation channel. This is due to 
the fact that the one photon annihilation part of the interaction has no additional 
singularity that needs to be taken care of numerically, because of the simple energy 
denominator ( [4 .21 ). We compiled our results in the form of binding coefficients in Table 
(|4.3|) . We have used there, contrary to previous chapters, the bracket convention for 
the errors to suppress the zeros. 

One notes a slightly larger breaking of rotational symmetry. The triplet ground 
states of different Jz are still approximately degenerate, but the discrepancy is bigger 
than without the annihilation channel. The dependence of these discrepancies between 
corresponding eigenvalues for J^=0 and Jz=^ on the number of integration points is 
shown in Fig. ([4.5| ). The behavior of the curves is similar to those of the calculations 
without annihilation channel. Fig. ( p.5| ). An additional plot. Fig. ( [4.5| ), shows the 
dependence of the rotational symmetry breaking on the cutoff A. The discrepancies of 
corresponding eigenvalues are almost independent of the cutoff A. 

To make a comparison of our results to those of perturbation theory easier, we 
show in Fig ( [4.7| ) the eigenvalues for a principal quantum number n=2 for both theo- 
ries graphically. The structure of the two plots is almost the same, only the 2^5*0 state 
and the 2^Pq state are exchanged in our results. This is due to the cutoff dependence of 
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AW r 




Figure 4.5: Deviation of corresponding eigenvalues for Jz — and = 1 including the annihilation 
channel (a — 0.3, A = 1.0m/): (a) as a function of the number of integration points N, (b) as a 
function of the cutoff A. The graphs show AM^ := = 0) - Af2( = l) for the states l^Si (A), 

2^ Pi (□), and 2^ Pi (o). 



the S-state, which is larger than that of the other states (cf Fig. ( |4.6| ) (b) and next para- 
graph). We stress that the results of the perturbative calculations change considerably, 
when the next higher order in a is considered. For example, the mass squared of the 
triplet 2^ Si is, according to Fulton and Martin |T|, M^{2^Si) = 3.9780186070 up 
to order a^Ry , which is nearer to our result than the value of perturbative calculations 
up to order 0{a'^) as displayed in Fig. (f4.7|). 



4.4 Parameter dependence of the spectrum 

The convergence of the eigenvalues with growing number of integration points is 
the same as the case of no annihilation channel. To be explicit, the convergence of 
the eigenvalues can be shown to be exponential. One can fit the singlet ground state 
eigenvalue excellently with the function 



M^{N) = M^{21) - M2(21)-M2(5) 



g-(iV-5)/2.2_ 



We did not perform the limit N ^ oo, because the accuracy of the results for > 20 
suffices to compare to other data. 

The cutoff dependence of the positronium spectrum including the annihilation chan- 
nel is comparable to that of the spectrum without it. However, a striking difference 
occurs: the inclusion of the annihilation channel stabilizes the cutoff dependence of the 
eigenstates. In particular, the triplet ground state in Fig. ( ^.6| )(a) shows only a small 
dependence on the cutoff, when one compares it to the behavior of the same state in 



Fig. (|2.8|) . One can fit these curves with a polynomial in log A. The singlet ground 
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Table 4.2: The binding coefficients of the singlet (Ss), the triplet 
(Bt), and the hyperfine coefficient C^j are listed for a = 0.3, 

Nx = 41,A^2 = 11- 



state eigenvalues are the same as without the annihilation graph and for the triplet one 
obtains 

^IrivieA^) - 3-91392 - 0.000288079 log A + 0.000147268 log^ A. 

Comparison with Eq. (|2.16| ) shows that the decrease of the triplet with log A is sup- 
pressed by including the annihilation channel by a factor of 60. Also the excited states, 
n=2, show a different behavior, as compared to Fig. ( p.8| ). Here, only one state shows 
level crossing as A grows large. The eigenvalue of the 2^Pi state, however, depends 
only weakly on the cutoff here. 

The values for the binding coefficients B^ and the coefficient of the hyperfine split- 
ting Chf are presented in Table ( [4.2| ). The values are correct for a cutoff A ^ 1m f 
when compared to results of perturbation theory up to (9(q;^). However, the effects of 
higher order correction to perturbative calculations are significant for a large coupling 
such as a = 0.3. The result of perturbation theory to 0[a^ In a) for the coefficient Chf 
is considerably smaller than that to order 

Concluding, we state that the cutoff dependence of the spectrum is improved as 
compared with the case of missing annihilation channel. This makes it evident that the 
annihilation channel is a necessary part of the theory. 

As a last investigation of the properties of our model, we will vary the coupling 
constant and interpret the spectrum. A similar procedure was performed by Dyks- 



HOORN ET AL. who studied coupled integral equations for QED-bound states in 
equal-time quantization with a variational ansatz. They calculate masses for the lowest 
eigenstates of positronium with and without the annihilation channel and plot them 
versus the coupling constant. The prominent feature of their figures is the occurrence of 
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Figure 4.6: The spectrum with annihilation channel as a function of the cutoff A: (a) ground 
states (ri=l), (b) first radial excited states (n=2). The parameters for the calculation are a — 0.3, 
Jz = OjNi = 2 5, N 2 — 21. One notices a better behavior of the curves with growing A compared to 
Figs. ( pTq ) and (2_^): the decrease of the eigenvalues is smaller, especially for the ground state triplet 
l^^i. Moreover, there are no level crossings for n=2, except the crossing of the singlet S-state at 
A = 1.5 (not visible in this plot). 



a critical coupling at which the masses become smaller than zero. We have performed 
the analogous calculations within our approach. It seems at first glance [Fig. ( [4.9|) (a)] 
as if the eigenvalues, after decreasing quadratically with the coupling as expected by 
perturbation theory, stabilize at a coupling a ~ 1.5. However, further investigation 
shows that this is merely an effect of the cutoff dependence of the spectrum. The 
eigenvalues in Fig. ( [4.9|) (a) were calculated with a cutoff A = 1.0 rrif, which is too 
small for an coupling constant of a = 1.0, since then the Bohr momentum is of the 
same order as the cutoff. Fig. ([4.9| )(b) shows clearly that there is a critical coupling. 
The masses calculated with a cutoff of A = 20 m/ tend to zero at a ~ 0.5. A similar 



value was found in |^ Ref. 2]. 



53 



Chapter 4. The annihilation channel 




CV2 



CO 
CO 



CO 



CD 
CO 



Figure 4.7: 
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Comparison of multiplets for n~2: (a) results of the present work with a 
A^2 = 21; (b) equal-time perturbation theory (ETPT) up to order 0{a^). 
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Figure 4.8: Compiled spectra of positronium with different Jz = —3, —2, . . . , +3 including the an- 
nihilation channel. All spectra have been calculated with a = 0.3, A = 1.0m/, Ni = N2 = 21. The 
mass squared eigenvalues in units of the electron mass m'j are shown. The notation for the states 

is ^^'^^Lj' . The resolution of the plot is inadequate for the multiplets. Nevertheless, the numerical 
degeneracy of the three triplet ground states ^S^^j^Sf, and ^Sl becomes very clear. 
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Figure 4.9: The spectrum as a function of tlie coupling constant: (a) A = 1.0 m/, (b) A = 20.0 m/. 
In (a), the eigenvalues seem to converge to a stable value as a grows large. However, this is just an 
effect of the cutoff dcipcndcnce of the spectrum: for a larger cutoff A = 20 m/, (b), the eigenvalues 
become negative at a critical coupling olc — 0.5. 
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Table 4.3: The positronium spectrum for a = 0.3, A = 1.0m/,iVi ^ N2 = 21. The non-relativistic 
notation for the terms and the eigenvalues for Jz=0 and Jz = + 1 including the annihilation channel are 
shown. The discrepancy between the eigenvalues is AB„ := Bn{Jz=0) — Bn{Jz= + !)■ The numerical 
errors arc estimated from the difference between the values for maximum and next to maximum 
number of integration points. The actual errors are smaller due to the exponential convergence of the 
eigenvalues with A''. The k numbers in brackets are the errors in the last k digits. 
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Chapter 5 



On the theory of effective 
interactions 



5.1 Introduction 

The aim of this chapter is to show that fixing the redundant parameter uj to the 
expression ( |2.1CID is not an ad hoc assumption, but a consequence of the structure of 
the theory and of the projection method used. The material presented here follows the 
line of arguments given in PO], where the method of iterated resolvents was introduced 



and applied to QCD. The case of QED is discussed in this section. 

To understand the following, it is helpful to review the reasoning presented in Chap- 
ter 1^ and to see things in a more general way. In the description of the positronium 
model it was pointed out that the unrestricted Fock space can be divided arbitrarily 
into two parts, one called the P- the other the Q-space. Since only two Fock sectors 
(|ee), 1 667)) were used in the model, it was clear how this separation had to be done. 



Tamm ||6^ and Dancoff |Q used this method in a different context in the follow- 
ing way. First, they truncated the Fock space to two sectors. Second, they projected 
one Fock sector onto the other. Third, the emerging energy denominator is modified 
to simplify the calculations. The so-defined procedure fails completely in front form 
dynamics p4| if the third step is missing. A severe (coUinear) singularity occurs that 



cannot be treated within this approach. Tamm did not recognize this problem, because 
he considered the modification of the energy denominator rather as a simplification than 
as an approximation for the effects of the omitted higher Fock states. 

What was done in Eq. ( p.5| ) was to introduce a new, redundant parameter u for 
the a priori unknown mass (squared) eigenvalue M^. This parameter is free and can 
be fixed on an ad hoc basis to remove the collinear singularity. It will be shown that 
this fixing, Eq. (|2.10| ), is not only the natural choice for this parameter, but that it is 
a consequence of the structure of the formalism. 
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5.2 The method of iterated resolvents 



The essential cause for the failure of a "naive" Tamm-Dancoff approach is the fact that 
the Fock space truncation is performed too early in the formalism. Because of this, 
one throws away all possible interactions with the omitted sectors. This gives rise to 
the problems described above. It is therefore necessary to study the structure of the 
resolvents to all orders and to investigate possible approximation schemes consistent 
with the solubility of the whole problem. 

Of course, one has to construct a finite dimensional Hamiltonian from the La- 
grangian density in the style of the Tamm-Dancoff approach. But a truncation cannot 
be the first step. Rather one introduces unphysical parameters needed to map the 
Hamiltonian operator onto a matrix in such a way that they remain in the formalism 
until the end. One can then investigate rigorously the behavior of the theory in the 
limit when all unphysical parameters are removed. 

We consider in this chapter the Hamiltonian within the DLCQ formalism, i.e. with 
discretized momenta. This maps the Hamiltonian operator into a matrix with a denu- 
merable number of columns and rows. It may still be an infinite dimensional matrix. In 
a Fock basis, the Hamiltonian is naturally decomposed into sectors of different particle 
content, cf. Table ( [4.11 ). Each sector by itself contains an infinite number of states, and 
is regulated by a cutoff A, like in Eq. ( p.l3|) . The number of Fock sectors, on the other 
hand, is limited by the finite longitudinal momentum P+, or rather by the harmonic 
resolution K = ■ The finite dimensional eigenvalue problem reads thus: 



( (0|if|0) (0|i/|l) 
(l|i/|0) {l\H\l) 

\ {N\H\Q) {N\H\l) 



m\N) \ 
m\N) 



{N\H\N) I 



\ km) I 



UJ 



( {m \ 



(5.1) 



Like in Chapter 0, one can reduce the dimension of the Hamiltonian matrix by project- 
ing the "highest" state {N\il)) onto the others. One arrives at an effective Hamiltonian, 
depending on the parameter a;, due to the resolvent 

Gn{uj) ■■= {n\ uj-H \n), 1 <n<N. 

In general, one has the eigenvalue problem 

n 

J2{t\H,X^)\j){j\ijiuj)) = E{u){tmu;)), 

j=0 

where the last state can be expressed by all others 



n-l 



(iV|^(u;)) = G^iuj)Y.{t\H^iu)\j){jmuj)). 

j=Q 
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Consequently, the effective Hamiltonian in the (n— l)-sector becomes for each matrix 
element 

= H^iuj) + H^{uj)G^{uj)Hn{uj), (5.2) 

which in a sense is a recursion relation. 

One has to find a systematic way to express the effective Hamiltonian after the z-th 
projection in terms of the bare Hamiltonian H. If one applies the relation Eq. (5.2) 
repeatedly, one obtains 

N 

Hn = H + ^ HmGmHm, (5.3) 

m,=n+l 

with the bare Hamiltonian H. For example, for n = 3 one gets 

H^ = H^ + H^iG^iHs + H2G2H2. (5.4) 

by inserting H2 = + H^G^H^ into Hi = H2 + H2G2H2. The general case can be 
proven by induction. An illustrative example is given in Section |5.4] . 

In Eq. (|5.3| ), the effective Hamiltonian in Fock sector n is expressed in terms of the 
bare Hamiltonian and scatterings into higher Fock sectors. It is important to notice 
that no scattering into lower Fock sectors occurs. From Eq. (|5.4| ) one can infer how 
the general structure of the expression for the effective Hamiltonian will be. Chains 
of terms with a different number of resolvents will emerge. In fact, it turns out to be 
much better to classify those chains by the number of its resolvents than, for example, 
by the order of the coupling constant. If all chains with k resolvents are collected in 
H'^'^\ the effective Hamiltonian is the sum 

i/„ = /fW + ^(i)+^(2) + ... . (5.5) 

This expansion is finite, because the bare Hamiltonian is a finite matrix. With this 
classification scheme, one obtains a recursion relation for the general term by inserting 
of Eq. (IH) into Eq. (|J) 

l>n 

The first terms read 



(5.6) 



l>n 
l>n 

Note the change in the order of the calculation. At first, the rows of the Hamiltonian 
were projected from to 1 unto one another to derive the expression for the effective 
Hamiltonian in the lowest sector. Now, the calculations are reversed in the sense that 
the effective Hamiltonian in the lowest sector is evaluated by longer and longer chains 
of resolvents, until one reaches the dimension of the bare Hamiltonian. The advantage 
is obvious: one is able to control the limit in which no truncation is made at all, 
i.e. N ^ CX3. 
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Figure 5.1: The three effective graphs of gauge theories (Courtesy of H.-C. Pauli). 



5.3 QED treated with iterated resolvents 

So far, the general method of how to calculate an effective Hamiltonian out of a fi- 
nite dimensional Hamiltonian matrix was described. In particular, no assumptions on 
the matrix elements of the starting Hamiltonian (like the one in Eq. ||5.1|| ) were made. 
From the structure of the QED Hamiltonian, described in Appendix |C| and displayed 



in Table ( |4.1| ), it is obvious that most of the matrix elements will be zero, certainly in- 
cluding all those in which the number of partons created or destroyed by the interaction 
exceeds two. This can clearly lead to significant simplifications. The procedure can be 
simplified even more by applying a technique of light-cone perturbation theory. There, 
one takes care of the instantaneous gauge parts of the Hamiltonian, i.e. of the seag- 
ulls and forks, only at the end of each calculation by redefining the non- instantaneous 
propagators Appx. A, Fig. 31]. One adds to the propagators a part containing an 
instantaneous graph including a ^-function for the longitudinal momentum transferred. 
One can apply the same trick here: one only keeps the vertex interaction in the matrix. 
Table (^). 

Since we want to calculate mesonic, i.e. fermion-antifermion systems, the sector |ee) 
plays the role of a "cornerstone" sector. If one calculates the effective Hamiltonian in 
this sector one gets 



Hgg = Hi = Ti + VGsV + VGiVG2VGsV, 



(5.7) 



where Ti is the kinetic energy in the |ee)-sector. This is an expression with at most 
three resolvents Gj{uj). The graphs of these fundamental chains are shown in Figure 
( |5.1| ): the photon exchange graph, the two-photon annihilation graph, and the self 
energy graph. 

It is worth mentioning two things here. Firstly, the resolvents are in general non- 
diagonal operators. Secondly, a single chain can symbolize more than on graph. For 
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example, the photon exchange graph and the self energy graph emerge out of the same 
chain VG^V. 

The structure of the effective interaction for the |ee) states seems to be very simple. 
The question arises what the effect of the higher Fock states will be. A first hint can 
be gained by looking at sectors similar to |ee), i.e. to |ee7), |ee77) . . . , 

: T^ + VGqV + VG(sVG^VGqV + VGiV, (5.8) 

- % + VGiqV + VGwVG^VGiqV + VGjV, (5.9) 

- T^^ + VGi^V + VGi^VGiiVGi^V + VGiiV. (5.10) 







Heeyf 




-^66777 - 


= Hio 



The only difference of these expression to Eq. (|5.?| ) is the absence of the last term, the 
annihilation of a photon into an electron-positron pair. This is because the first sector 
does not contain a photon. 

These sectors seem to be similar to each other. However, they are only a small set 
in the whole Hamiltonian. Their defining property is that the photons do not interact 
with the fermion-antifermion system. To classify the sectors, this kind of interactions 
is named spectator interaction Un- The other class of interactions, where the photons 
do interact, is called participant interaction Un. Consequently, the decomposition of 
the Hamiltonians read 

Hn = Tn + Un + Un , n = 3, 6, 10, 15, . . . . (5.11) 



Analogously to Eq. (|2.5|) , where the resolvent was expanded around its diagonal 
{i.e. free) part, the resolvents Gn can be expanded around the part containing the 
spectator interaction [/„ 

Gn = 7^-^- (5-12) 

UJ - In - Un 

The full resolvent is obtained by the infinite series 

Gn = Gn + Gn Un Gn = Gn + GnUn Gn + GnUn GnUn Gn + . . . . (5.13) 

Note that the unperturbed resolvent contains the interaction! Moreover, one observes 
that in this approach, the system does not leave the sector n. 

Let's repeat what we have achieved up to now. We divided the interaction in one 
part that acts only between the electron and the positron (spectator interaction) and 
one part where the photons can interact with the two fermions (participant interaction). 
We expanded the full resolvent Gn as an infinite series around the spectator interac- 



tion. One can show pO|, that the spectator interaction allows for three fundamental 
graphs only. Fig. ( ^.1| ), if one reduces the bare Hamiltonian by subsequent projections 
to the fermion-antifermion sector. This is the end of the procedure for non-abelian 
gauge theories. In QED, however, one has an additional sector I7) responsible for the 
annihilation diagram. Fig (|4.2| ). The whole spectator interaction resides in the coupling 
function attached to the vertices. 

The question arises, if the spectator interactions in the different sectors are somehow 
related. Note that the only difference between their diagrams is the different number 
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of non-interacting photons. In fact, this is the crucial point in the formalism. We will 
give here a heuristic explanation, how this can be understood. 

Consider the separation of the Fock space into two parts, named P- and Q-space, 
like in Chapter According to the discussion of this chapter, the only interaction in 
both spaces is the spectator interaction. What does the resolvent look like? It is from 
Eq. (U) 

Giu) := , (5.14) 

where Hq = Tq + Vq is the sector-Hamiltonian of the Q-space. The redundant param- 
eter u stands for the actual eigenvalue of the whole eigenvalue problem and contains 
therefore a kinetic and an interaction part 

^ = ^true + ^rue- 

The resolvent reads thus 

G{u;) = ^ -. (5.15) 

true + I'true — J-Q — VQ 

With the help of two smallness assumptions on the momenta of the photon 

< 1, and < M^g. (5.16) 



it has been shown in 20 , that indeed 



lim Vtrue -Vq = Q, 

which is very plausible from what we have stated before: the interactions in the different 
sectors deviate only by a different number of non-interacting photons. We stress the 
importance of the continuum limit at this point. If one has a finite number of Fock 
sectors, the argument does not hold. We obtain finally 

true ~ -'■Q ^ 

By this, one has shown that the resolvents are diagonal in the solution if one makes 
the two smallness assumptions ( ^.16|) , which are fully justified within a bound-state 



formalism. In fact, one makes the surprising observation that the resolvents are totally 
independent of the redundant parameter uo. It is clear that this parameter is merely a 
mathematical tool to perform the calculations in a controlled way. 

We recall the result of the plausibility argument in Chapter |^, namely the definition 
of the energy denominator P, Eq. (|2.9|) , together with the analytic expression for T*, 
Eq. ( |2.1CI| ), in the |ee7)-sector. Evaluating Eq. ( |5.17| ) in the |ee7)- and |ee)-sector, one 
gets by the arguments of the general formalism exactly the expression of Chapter |^, 
Eq. dnO). 
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One might wonder why the correct fixing of the redundant parameter a; resuhing in 
the diagonal resolvent Eq. ( |5.17|) could appear as an approximation of the expansion 



Eq. ( p.8|) , or even as an ad hoc assumption in |23, Eq. (2.7)]. The answer is that it is 
very important for the theory, at which the stage the (necessary) truncations are made. 
If one truncates too early in the formalism, it is clear from the work of Tamm and 
Dancoff that other prescriptions, such as ad hoc assumptions, have to guarantee 
the solubility of the equations. This is handled better when the truncation happens in 
a controlled way as in the formalism of iterated resolvents |^ . 



The general formalism of effective interactions can be understood as a summation 
over all intermediate states in the effective fermion-antifermion sector. It is, however, 
no Tamm-Dancoff truncation, as sometimes referred to in the literature. In fact, the 
Coulomb potential comes out of the formalism correctly only if one sums over all 
photons and fermion-antifermion pairs, i.e in the limit iV — oo. 

Note that the renormalization problem remains unsolved. After having derived the 
effective interaction within the present formalism, the cutoff dependence of the results 
have to be investigated. This is still a "stumbling stone" for the case of QCD. In QED 
the problems are present, but seem to be not as dramatic. 



5.4 Application to a model Hamiltonian 

To give an instructive example of how the method of iterated resolvents works, we solve 
the eigenvalue equation 

H^°y\'4>) = E\i>) 

where the 5x5 dimensional toy Hamiltonian is 
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(5.18) 



with D := 10, ^ := 4, 5* := 3, F := 1. The eigenvalues Ei can be calculated by standard 
methods 



Ei E {4.1386,6.9387,9.4239,11.0368,18.4908}. 



(5.19) 



The space in which this matrix operates is divided into two subspaces |1) and |2), 
i.e. N=2 in Eq. ( [5.1| ), so that the diagonal block matrices are 
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The effective Hamiltonian in tlie |l)-sector, according to tlie method of iterated resol- 
vents (cf. Eq. ||5.3|| ), is 



N 



ii>i 

where the resolvent 
Go 



1 



u - {2\Ht°y\2) 



j u-D V 
V u- D 
V u-D 



is itself a 3 X 3 matrix. The Hamiltonian to diagonalize is a 1 x 1 block matrix with 
two columns and rows. The corresponding two eigenvalues, both functions of cu, are 
plotted in Fig. ( [5.21 ). Fixing the redundant parameter uo by the condition 

E{uj) = u (5.20) 



yields exactly the eigenvalues given in Eq. (|5.19|) . Note that the eigenvalue functions 



for the first and second eigenvalue have continuous transitions. This is due to the 
diagonalization algorithm used which will always consider the lowest eigenvalue to be 
the first one. Note also the strange behavior of the curves yielding the second and 
third eigenvalue of the whole matrix. One sees that the eigenvalues repel each other 
in the numerical program. However, it is clear that there should be a crossing of the 
two curves. In fact, we have two functions in the plot, one of which has one and the 
other has two poles. This yields five intersection points with E{ijj) = u, as expected. 
Another important observation is the fact that the dependence of the functions on u 
is rather weak at the intersection points, which yields numerically stable eigenvalues. 

The result obtained should be compared with the similar plot in Fig. 5] . There, 
a 4 X 4 matrix was solved by applying the projection method row by row, and one 
ended up with an effective 1x1 Hamiltonian, a function of u, but nonetheless a real 
number. Here, we have shown that the formalism can applied also to matrices: the 
effective Hamiltonian is a matrix depending on uo. But the functions for the different 
eigenvalues Eiiuj) combine in such a way that continuous functions emerge which have 
N intersection points with the line defined by Eq. ( |5.2U| ). This corresponds to the N 
eigenvalues of the toy matrix ( |5.18| ), just as in the case discussed in [pT| . 
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Figure 5.2: Eigenvalues of the effective 2x2 dimensional toy Hamiltonian as functions of the 
redundant parameter uj: (a) the first, (b) the second eigenvalue, (c) the two plots in one figure. Note 
the continuous transition from one eigenvalue to the other. The five eigenvalues at the intersection 
points with the fixing condition E{uj) = uj (dashed line) are encircled. 
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Chapter 6 

Summary and Outlook 



In this thesis, the calculation of the spectra and the wavefunctions of an electro- 
magnetically bound, relativistic fermion-antifermion system in (3 + 1) dimensions has 
been performed. The problem has been solved for all components of the total angular 
momentum, J^, and for an arbitrarily large coupling constant. The latter required the 
application of non-perturbative methods in quantum field theory. An effective Hamil- 
tonian in the |ee)-sector has been used. The corresponding integral equation is mapped 
into a matrix eigenvalue problem and solved numerically. The spectrum and the wave- 
functions of positronium have been calculated in the regime of an unphysically strong 
coupling constant, a = 0.3. A computer code has been created that yielded numerically 
stable results with a minimized demand on computer equipment. The numerical meth- 
ods described in ["M] have been tested and improved. In particular, the counterterm 
technology for the Coulomb singularity has been fully implemented. Since the coun- 
terterms may in general involve an integral which is not analytically solvable, one can 
use only numerical integrations and the numerical effort increases. As can be seen by 
comparing the results of this work with p3], the convergence is improved noticeably. 



The most prominent aspect of the numerical results is that the spectrum is found to 
coincide to a very high degree of accuracy with the results of elaborate perturbation 
techniques. Previous work in this direction (Tang [^, Kaluza [^) failed to produce 
significant results because of numerical problems. 

Special care has been taken in the comparison of the eigenvalues of the positronium 
model used with results of perturbation theory. In the sector Jz=0 the Hamiltonian was 
split up into four block matrices with definite quantum numbers of charge conjugation 
and T-parity. This allows a one-to-one comparison of states with those resulting from 
perturbation theory, since the symmetry quantum numbers can be related to those 
used in instant form dynamics. In addition to the fact that the multiplets of a fixed 
principal quantum number n contain the correct number of states up to at least n = 5, 
we find that the correct combinations of quantum numbers are included in each of these 
multiplets. A careful study of the higher excited states shows that the calculations of 
this work are trustworthy even for these states, as opposed to, for example, lattice 
gauge theories. 
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Chapter 6. Summary and Outlook 



As was pointed out, rotational invariance is problematic when using front form 
dynamics. The operators attached to the transverse rotations are complicated and 
contain the interaction explicitly. It is therefore important to study the full spectrum 
of positronium by generalizing the model to a non- vanishing component of the total 
angular momentum. We find that the matrix elements can be calculated even for arbi- 
trary Jz- The improvement of the counterterm technology was crucial in this case. It is 
important that the counterterms are taken care of individually, because the singularity 
structures of the diagonal matrix elements differ considerably from each other. We 
stress the fact that this counterterms are due to numerical problems and have nothing 
to do with the physics of the positronium model. We neither use any special prescrip- 



tion to restore gauge invariance, as suggested by Coester [TG], nor do we construct 



non-covariant counterterms for this purpose, as necessary in front form perturbation 
theory |5l||. One way to study rotational symmetry in front form dynamics requires 
the diagonalization of operators at least as complicated as the light-cone Hamiltonian 
itself. We listed these operators in the first section of Chapter ^. It is clear that this is a 
non-trivial task. Almost all efforts up to date are spend to construct and to diagonalize 
the light-cone Hamiltonian, and not the other dynamical operators. Surely, the latter 
is an interesting project. However, we were able to construct an effective Hamiltonian 
and to solve for its spectrum in all sectors of different z-components of the total angu- 
lar momentum. In this way, we can learn from the physical observables (mass squared 
eigenvalues), how quantum numbers concerning the complicated operators can be at- 
tributed to the states. We point out the importance of the non-perturbative approach 
used in the present work. It has been shown in several models (even QED3+1) that 
perturbation theory in front form dynamics yields results where rotational symmetry 
is explicitly broken. There are no degeneracies of states with the same total angular 
momentum. This is clear from our point of view. Exact rotational symmetry is ex- 
pected only if it is summed over all Fock states and all momenta. Because we used an 
effective theory in the present work, where the effects of all Fock states are included in 
the effective interaction, we observe the desired degeneracies in our approach. 

The properties of the spectra calculated for different values of are the following. 
Most remarkably, we find that states of different but same total angular momentum 
J are numerically degenerate. There are slight deviations from exact degeneracy also 
for a large number of integration points. The latter can be explained by the restrictions 
of the positronium model used, e.g. finite cutoff and others. These degeneracies make 
it possible to classify the states according to their quantum numbers of total angular 
momentum a posteriori and make the diagonalization of the complicated operators of 
transverse rotations superfluous. Even more, it shows that rotations are unproblematic 
on the light cone, because rotational symmetry is restored in the solution, a previously 
worrying fact. The results show very good agreement with the values obtained by 
equal-time perturbation theory up to order 0{a^). 

The model was enlarged by including the one photon annihilation channel. The 
inclusion is a test for the consistency of the model. In particular, the conception of 
an effective theory operating with resolvents can be falsified if the annihilation channel 
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cannot be implemented without special assumptions. Our results show that the im- 
plementation of the annihilation channel is unproblematic in the sense that the same 
formalism can be applied as in the case of the projection of the effective |ee7)-sector 
onto the electron-positron sector. As an interesting property of the annihilation channel 
we find a strict separation of the instantaneous and the non- instantaneous interaction: 
the seagull interaction is present only in the Jz = sector, whereas the dynamic graph 
has non-vanishing matrix elements only for = ±1. Our approach passes another test: 
both graphs yield the same contributions to the eigenvalues and consequently rotational 
symmetry is seen also in the spectrum including the annihilation channel. Moreover, 
the inclusion of the annihilation channel improves the results for the hyperfine splitting. 

We stress the point that the implementation of the annihilation channel completes 
the investigations of how to construct an effective interaction for the electromagnetic 
fermion-antifermion system in the meaning of the method of iterated resolvents, de- 
scribed in Chapter ^. We have put all effects of higher Fock states into an effective 
|ee7)-sector, as far as the spectator interaction is concerned, i.e. the interactions in 
which the photons are not directly involved. The remainder of the interaction relies 
in the coupling function of the vertices. A hint to this conclusion is the logarithmic 
cutoff dependence of our results. It is clear that the couphng constant depends on the 
cutoff and has to be analyzed. A future aim, beyond the scope of this work, will be to 
show that the physical results of our model become independent of the cutoff, as soon 
as renormalization group techniques are consistently applied. This is supported by the 
fact that we find a stabilizing effect of the annihilation channel on the dependence of 
the spectrum on the cutoff: all eigenvalues show slower variation with growing cutoff 
A when the annihilation graph is added, in some cases this even prevents level cross- 
ings with other states. We therefore conclude that our model is correct as long as the 
vacuum polarization effects are not considered. 

The possibility of one boson annihilation is the main difference between QED and 
QCD in effectively truncated Fock spaces where the dynamic three- or four-gluon in- 
teraction is not possible but resides in the coupling function. The one boson sector 
exists only in QED because of the constraint of color neutrality in QCD. We showed 
that this sector can be projected onto the (already) effective electron-positron sector. 

In Chapter |^ we described the method of iterated resolvents which allows for the 
construction of effective interactions on the light-cone. We applied this formalism to 
our problem. It was shown that the assumption is correct that the effects of the higher 
Fock states can be accounted for by the fixing the redundant parameter by a function 
of the light-cone momenta, Eq. ( p.lO ). The fixing was necessary to derive an effective 
integral equation starting from the QED(3+i) Hamiltonian and is equivalent to the 
calculation of the effective resolvent in the |ee7)-sector. We showed that this resolvent 
is indeed diagonal in the Fock basis used and is therefore a function rather than an 
operator. This fixing procedure was supported before by the numerical precision of the 
calculations using it as an assumption. 

Because of a confusion of technical terms we stress the following in this context. 
The formalism applied here is misnamed when referred to as a Tamm-Dancoff approach. 
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Chapter 6. Summary and Outlook 



The Tamm-Dancoff approach was invented for nuclear forces in instant form dynamics. 
Its three main steps become very clear in the work of Tamm [|67| , §2]. The first step is 
to truncate the Fock space to two sectors (P- and Q-space). Then a projection of one 
space onto the other is performed. The emerging energy denominator (or resolvent) 
is modified to render the system solvable as a third step. What we have done is (a) 
regulation of the number of Fock sectors, (b) derivation of an effective Hamiltonian 
by subsequent projections of Fock sectors onto another, (c) proof that in the limit of 
infinitely many Fock sectors, the resolvents become diagonal for bound states. In our 
approach it is clear that a sheer truncation must fail, because essential parts of the 
Fock space are missing. It remains to treat the so-derived effective Hamiltonian with a 
renormalization group analysis. We have shown that the dependence of our results is 
logarithmic (weak) in the cutoff A. The inclusion of the annihilation channel stabilizes 
the eigenvalues. Whether or not this hints to a weaker dependence of the (running) 
coupling of QED on the cutoff when compared with the coupling of QCD, where this 
channel is absent, is subject to further investigation. 

This thesis was motivated by the will to understand the problems of non-perturbative 
methods on the light-cone. The hope is clearly to find an adequate description of the 
low energy region of QCD and thereby of the QCD bound states, like mesons and 
hadrons. Front form dynamics offers very promising features, such as a simple vac- 
uum, and more. Other methods have been proposed to deal with these problems, 
namely approaches which diverge from the constituent quark picture. The method 
of similarity transforms of the bare Hamiltonian proposed by Wilson et al. 
and the coupling coherence technique advocated by Perry et al. seem to be very 
promising. Although the formalism is very elaborate, there are problems in calculating 
actual numbers. A model of Brisudova et al. |^ aims at calculating mesons with 
one heavy constituent quark with these methods. However, the authors are forced to 
perform a line of approximations, and end up with a non-relativistic, rotationally non- 
invariant Schrodinger equation, which makes it hard to tell the actual achievements of 
the approach apart from the draw-backs of the approximations made. An instructive 
analytical study is that of Jones et al. |8^, which makes the approach of Wilson 
and Perry very clear. The authors calculate the positronium spectrum up to order 
0{a'^) and show that the results are independent of the cutoff up to order 0{a^). The 
triplet ground states are degenerate. This is not surprising because the effective inter- 
actions are derived for small couplings a -C 1 and a non-relativistic limit of the theory 
is considered. It is known (Eq. p.l9|) that the term that breaks rotational invariance 



vanishes in the non-relativistic limit. A similar method is that of Wegner 84 , who 



uses flow equations to renormalize the bare Hamiltonian. Attempts have been made 



B5| to apply this formalism to QED and to calculate the positronium spectrum. A 
link between the DLCQ approach and the method of Wilson et al. is the work 
of Ammons There, a Tamm-Dancoff truncation is performed together with an 

application of the similarity transforms. These approaches are distinct from the one 
used in the present work in the following sense. First, the Hamiltonian is constructed 
using renormalization theory, then the associated eigenvalue problem is solved. Up 



70 



to date, the latter step has been performed using bound state perturbation theory or 
other perturbative methods, ahhough in principle a non-perturbative solution seems 
possible. Complementary to this formalism, we construct first an effective Hamiltonian 
to account for the many-body aspects of the theory |8^, solve for the spectra with a 
non-perturbative method and then have to analyze the results using renormalization 
techniques to study the implications of Quantum Field Theory. 

Concluding, the formalism of DLCQ together with the theory of effective interac- 
tions, as described in Chapter ^, can be combined to yield non-perturbative results 
for a relativistic bound system in very good agreement with other standard methods. 
Even in (3+1) dimensions it is possible to construct a consistent and solvable model 
for bound states with an arbitrarily large coupling. 

With the results of this work and the computer code at hand, one can now proceed 
to attack QCD bound states in true (3 + 1) dimensions. Because we have reduced the 
Fock space by effective methods, the difference between QED and QCD relies in the 
coupling constant, or rather coupling function. The method of iterated resolvents [pO[] 
allows for the calculation of these coupling functions for both QED and QCD, but this 
is tedious work and has not been tackled. On the other hand, Merkel et al. 
used a phenomenological running coupling for QCD, inspired by the work of Richard- 
son [0, to construct an integral equation, analogous to the one in the present work, 
for the bound states of a quark-antiquark system. This equation is solved with varia- 
tional methods. The fitted meson masses show good agreement with experiment. This 
phenomenological coupling can now be plugged into the formalism and computer code 
derived in the present work to solve for the spectrum of a QCD-bound quark-antiquark 
system. This seems possible, because the full effective Hamiltonian of a gauge theory 
contains only three essential graphs, which themselves include a coupling function. The 
computer code created in the present work already contains two of these graphs, and the 
third one (the two gluon annihilation) is known to be less important. Consequently, 
the first task is to find a reasonable phenomenological running coupling constant to 
construct such a test. The main problem is to find an appropriate regularization of 
the severe (p~^) singularities occuring. Because of this, the dependence of the results 
on the regularization scheme used has to be investigated carefully. Nevertheless, this 
seems a promising way to proceed to understand hadrons as QCD-bound states. 
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Appendix A 
Notation 



The notations used in this work are compiled in the sequeL New technical terms are 
italicized when first introduced. The conventions used for discretizing the theory can 
be found in Appendix ^. 

Coordinates 

The light-cone coordinates are defined by 



The greek indices run like (+,1,2,—), latin indices denote transversal directions, e.g. 
i = 1,2. Analogously the momentum coordinates of a particle are {p'^ , 'p±, p~) . Hence, 
p~ is the light-cone energy of a particle. An underlined variable denotes a merely spatial 
vector X := {x~,x±). Some frequently used symbols are: 



X 




We use the metric 



9' 



/ 2 \ 
0-100 

0-10 

V 2 / 



P+ 



total longitudinal momentum 

total transversal momentum 

longitudinal box length 

transversal box length 

harmonic resolution (integer valued) 

discretization volume 

summation over all quantum numbers 



2L 
2L± 



K := 



TT 



n := 2L{2L±) 
all QN 



The relative coordinates of the i-th particle are 




and 



k±i '.— XiP P±i} 



(A.l) 
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referred to as the longitudinal momentum fraction and the (relative) transversal mo- 
mentum, respectively. If Pi := 0, they have have the properties 

^Xi = l, and ^ki = 0. 



Commutation relations 



The commutation relations according to the Dirac-Bergmann (pQl [90] algorithm are 



A\x),d+A^{y) 



To be consistent with these relations, in the expansions of the fields the operator-valued 
coefficients have to obey 



All other (anti-) commutators vanish. 
Coupling constant 

We set 



a 



9_ 
Att' 



In QED: g = |eo|, a = 1/137.0359895(61) = 0.007297353. In the discretized theory we 
use 



P 



Spinors and polarization vectors 



The fermion fields are separated into two different helicity eigenstates by 
Here, the projection operators read 



1 



A± := -7 7 



or explicitly 



A. 



A_ 



(I -1 
10 1 
-10 10 

V 1 1 y 

The explicit form of the spinors and polarization vectors used is given in Eqs. (p3.5| ) and ( [B.6|) . 
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Appendix B 

QED(3+1) on the light-cone 



The Hamiltonian operator of Quantum Electrodynamics in (3 + 1) dimensions in front 
form dynamics is derived. QED describes the interaction of (electrically) charged parti- 
cles. The interaction is intermediated by the gauge field of the photons. We begin with 
the Lagrangian density of QED(3_|_i) 

ip — rufipip — gijj^^ipA^. 
The field strength tensor is 

The equations of motion emerging for the gauge fields are the Maxwell equations 

d,F>'^^f:^gi,Yi^, (B.l) 
and for the Fermi fields the Dirac equation 

{iYD^-mf)iP^O, (B.2) 

where 

Df,:= + igAf, 

denotes the covariant derivative. We use in this work the light-cone coordinates 

:= x° ± x^. 

x~^ plays the role of a time, x~ is a direction of space. To write scalar products involving 
Dirac matrices likewise, we define 

We use the usual Dirac representation for these matrices. We work in the so-called 
light-cone gauge 

A+ ^A° + A^ = 0, (B.3) 



r — F pfj-i" 



^7^9^ - (9^^)7^ 
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which leads to a consistent theory in the normal mode sector. Deriving the light-cone 
energy P~, we shall follow the line of arguments of Tang ||4^, Ref. 1]. With help of 
the projection operators 

A± := -7V = ^7^7^, 
we separate the fields into states of different hehcity 

^± = A±^. 

It turns out that only six out of the twelve fields A^^, ip^, and occuring 

in the Lagrangian density are independent of the others: the conjugate momenta of 
A~ , ip-, and ipl_ vanish. The Euler-Lagrange equations in these fields contain no time 
derivative, and have the rank of constraints. This is a consequence of the chosen 
quantization plane. The equations of motion for the dependent fields are 

id^Tp- = {^—idio' + [3mf + gAia^^ 
id^ipl = -ipl(^idia' + Pmf + gAia^, 
{id+^A- = 2d+diA^ + Ag^l'ilj+. (B.4) 

with the usual Dirac matrices a\(3 and di:= diipl^. 

To quantize correctly, we have to eliminate all dependent degrees of freedom. This 
can be achieved either^ by inverting the equations ( p.4|) or by applying the Dirac- 
Bergmann algorithm |^9|, |91], The latter leads for a reasonable Lagrange density 
necessarily to the right commutation relations for the independent fields. 

We substitute the dependent fields by functionals of the dynamic degrees of freedom 

= -^(-idia' + (3m})i)+- -^A'a'ilj+, 
= ^i'i{id,a' + l3mf^ + ^^\A'a\ 

The inverse derivatives are defined by Green functions^. From the energy-momentum 
tensor 

one obtains the momenta by integrating over all space directions 

P'^ = - / dx- £xs_T+''. 

2 J —00 J —00 

^Cf. |44[ Chapter 4]. 

^These definitions are not unique |^. 
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Appendix B. QED(3+1) on the light-cone 



The light-cone energy P has the following structure 



P' = - dx~ (fx^V^ + gV^ + g^V^. 

^ J -~oo J —oo 



The different terms are 



sym 



-2 V- 



-4 ^^' 



i di + (3m 
1 



id' 



sym 



{id+f 



id+ 
id^id,A^\ 



A'a^ip^ 



sym 



2 \ ^plAW—A^a^ij+ 



{id+) 



sym 



Here, the symmetric brackets are defined by 



sym 



A- 



^(9^ 



:5 



■.= A 



sym 



{id+y 



_ 1 

" 2 
r5 + 



Z(7+ 



-A 



id^ J \id^ 



1 



A\B 



:A B. 



We expand the system in plane waves at light-cone time = 0, and work subse- 
quently in momentum space. The notations used here are those of [Q. By restricting 
the system to a box 

-L<x'< L, 

we discretize the momenta. We have to impose boundary conditions compatible with 
the equations of motion ( p.l|JB.2D . The expansions of the dynamical fields with the 
notation x := {x~,x±) for space and k := {k^,k±) for momentum variables read 
explicitly 



ip+{x 



A\x) 



V ' ' s n 



^ A p vk^ 
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Here, Q = 2L{2L±)'^ symbolizes the volume of the box, in which the momentum vectors 
have to lie. 

We are forced to demand periodic boundary conditions in all space directions for 
the gauge fields, since they couple to a bilinear term. The boundary conditions for 
the fermion fields are not subject to such a constraint. It is convenient to impose 
antiperiodic boundary conditions merely for the longitudinal fermion field. The zero 
mode of this field vanishes by that choice. Of course, this is of no importance to us, 
since we consider a priori only the normal mode sector by working in the light-cone 
gauge. The summations are performed in the following way. For fermions the indices 
are 



and photons have 



nil 

T 
T7 



pir 



n = 1, 3, 5, . . . 

n' = 0,±1,±2, 

p = 2,4,6,... 

p* = 0,±l,±2,. 



The projected spinors are defined according to Lepage and Brodsky p3| to be 
M+(A) := and f+(A) := x(— A). The spinors 



x(T) = ^ 



1 

V2 



/1\ 


1 

V J 



1 

^72 



/ \ 

1 



V -1 J 



(B.5) 



obey the relations 



xnA)x(A') 



EXa(A)x^(A) = A, 



a/3- 



The polarization vectors 




and e(|) : = 



are orthonormal in the transverse space and complete 

e*(A)e(A') = 5a,a', 
Eei(A)e*(A) = 6ij. 



(B.6) 
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Appendix B. QED(3+1) on the light-cone 



The projected spinor ip+{x) can easily be generalized to a spinor for the full free fermion 
field ip{x) by substituting m+(A) and t'+(A) by 

u{k, A) = + f3mf + d±- x(A) 

vik,\) = ^= + /3m/ + a_L • x(-A). 



The same is true for the gauge field, if the purely transversal vector ej_(A) is substituted 
by the four vector 

/ 



k+ 

V 6-l(A) 



We quantize by postulating the usual commutation relations for the coefficients of the 
Fourier decomposition of the fields. The coefficients become operator valued and obey 

{b, b} = {d, d} = {b, d^} = [a, b] = [a, d] = [a, b^] = [a, d^] = 0. 

Since we address ourselves in this work to solve an integral equation, we additionally 
list the expansions of the fields in the continuum 



1 



1 



E 

A 

E 



dk^ 

dk^ 



oo ^ 

d'^k± — 

-oo y/k+ 

oc ^ 

d^k± — 
-oo Wk~^ 



b{k, \)u+{\)e"^-^ + d\k, \)u+{\)e+'^-^ 
a{k, \y{\)e-'^-^ + a\k, A)e«(A)e+*^-^ 
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Appendix C 



Matrix elements of the light-cone 
Hamiltonian in QED 

The matrix elements of the QED (3+1) Hamiltonian in light-cone quantization have been 
derived elsewhere [^. The light-cone Hamiltonian is defined as 

It is usually divided into the parts 

Hi^c = T + V + S + F + C, 

where the kinetic energy is 

r = E ^^^^ ("5", + + E v'^'H- 

q -l^q q -^q 

V, S, F and C are the vertex-, the seagull- and the /orA;-interaction, as well as the 
contractions [|. The Fork interaction is absent in this work because of the truncation of 
the Fock space and is not listed. 

The creation operators and create plane waves (states) for electrons, positrons 
and photons, respectively. These particles are characterized by their quantum numbers 

q := {x,k±,X). 

— * 

A particle with longitudinal momentum fraction x, transversal momentum k± and he- 
licity A is annihilated by the (annihilation) operators bq,dq,aq, depending on the nature 
of the particle, and represented in the graphs of the matrix elements by the following 
symbols 

^ This nomenclature stems from the operator structure of the matrix elements. Is the parton number 
changed by an interaction by 0, 1 or 2, the operators are called seagull, vertex, or fork, respectively. 
Contractions occur due to normal ordering of the seagull operators. 
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Appendix C. Matrix elements of the light-cone Hamiltonian in QED 



► = {x, k±, X) (electron), 

4 = d^{x,k^, X) (positron), 

^\/\/\/\/ aj = a^{x, k±, X) (photon). 

The photons are massless, the fermion mass is ruf. Creation and annihilation op- 
erators obey the usual commutation relations. It is summed over all possible quantum 
numbers in all expressions. 

The matrix elements are listed in the discretized form where the normalization 
volume is 

with the half longitudinal and transversal box lengths L and L^, respectively. The 
coupling constant in the tables is g"^ = S'^prn, the total longitudinal momentum P+ = 
j^K, with the harmonic resolution K. For all other notations see Appendix ^. The 
continuum limit is obtained by substituting simultaneously 




Graph 


Matrix Element 


' " !) 


• 1 


oo 


{ ' - ' ) 

\{xi — x'Y {xi+x'Y) 




1 — *— 


'J 


> 1 


oo / 


1 1 ^ 

^x'{xi + x') x'{xi—x')j 


1 






r\/^ 1 


oo 
x' k' 


( ^ 1 ^ 

\Xi{x' + Xi) Xi{x' — Xi] 


■) 


c 


all QN 


+ C<«>(l))+aWc<«'(l); 





Table C.l: Matrix elements of the contractions. For further explanations, see text. 
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Graph 



Matrix Element 



Helicity 



1 — 



K-..(i;2,3) 



irif / I 1 
+9^ 



+9\ —e±{X3) 

V 2^3 \X3 
[Y^ (k^3 

+g\ —e±{X3) — 

V 2^3 \ X3 



k±2 

— * 



X 




V-,_,,(l;2,3) = ~g 



X3) 



Xi \X2 



-^W— e± Ai) 

\ xi yxi 

-g\ 

\ Xi \ Xi 



k±3 
X3 

X2 



^ = E {b\b2a3-d\d2a3)Vg^gg{l;2,3) 



all QN 



+ E {alblh-aldld,)v;_^Jl;2,3) 



all QN 



+ E [4&24^;-.,,-(i; 2, 3) + 46^1 2, 3) 

all QN 



Table C.2: Matrix elements of the vertex interaction. For further explanations, see text. 
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Appendix C. Matrix elements of the light-cone Hamiltonian in QED 



Graph 


Matrix Element Helicity 


1 > f *— 3 

2 * * * 4 






SSI„-(1,2;3,4)=9^I^^-^^ xi\i\ 


2-" J 


/; 




1 » 

2 'S/V/N/V 


» 4 


S£L„(1.2;3.4) = ^^^_;^^^ 


^ = E ^I4M4[^£1,5(1,2;3,4) + ^J^-1,,-(1,2;3,4)] 

all QN 

+ (6l4M4 + 44c?3a4) [^^J_.,,(l,2;3,4) + ^(^L,,(l,2;3,4)] 

all QN 



Table C.3: Matrix elements of the seagull interaction used or mentioned in the present work. The 
full table can be found in B3]. For further explanations, see text. 
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Appendix D 

Calculation of effective matrix 
elements 



The calculation of the matrix elements of the effective interaction using matrix elements 
of the canonical Hamiltonian and using currents is described in short. We use the 
derivation of Krautgartner, Pauli and WOLZ [^. The interaction within the 
positronium system can be expressed via currents ||2^, Eq. (2.13)] 

V V \x — x'^ 



Here the currents are 



the momentum transfers read 

and the energy denominator is given by 

V =\x-x'\{T* -uj) -ll (D.2) 

using = i (^e + ^i)- order to derive ( p.l| ), one evaluates the diagrams shown in 
Fig. (p.l| ) according to the rules of light-cone perturbation theory as formulated by 
Lepage and Brodsky in Appx. A]. In the energy denominator, the average of the 
energies of the in- and outgoing particles is used instead of the sum of the energies of 
the incoming particles only. The second term of Eq. ( |D.1D vanishes, if one fixes u = T* 



as described in Chapter ^ The definition of the interaction matrix elements is read off 
from the integral equation. We follow ||3^, Eq. (A. 5)] and Eq. (3.19)] by setting 
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Appendix D. Calculation of effective matrix elements 



x{l — x) 



dx'd^k'^ {x, kr, Ai, X2\j{lerj{k)\x', A;l; s[, 4) 



xx'{l — x){l — x') 



1pn{x,k±;Xi,X2) = 0. 



For the actual calculation of the effective interaction via the currents we need the matrix 
elements of the Dirac spinors: 



M 



=u{k\X)Mu{k,X) 



r 



r 



k+k'- 



7i 



k'+ 



+ )^'' "^^[k'^ k+)-'' 



72 



k'- 



k+ 



k'+ k+'-'' 



Table D.l: Matrix elements of the Dirac spinors. 

With the notation 

u{x', k'j_, X[)Miu{x, k±, Ai)it(l - x', -k'j_, X2)M2u{l - x, -k±, X2) 



xx'{l — x){l — x') 



{M1M2) := 



the prescription for the calculation of the matrix elements of the effective interaction 
reads 



{x, kr, Ai, X2\f{le)Uk) \x', k'^; X[, A'2 



{x, kr, Ai, X2\j^{ie)jMW, fej; a;, a^^) 

2^xx'{l- x){l- x') 
^(^(7V) + ^(7V)-(7?)-(72'))- 



(D.3) 
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1-x.-k^ 1-x',-k'j 



+ 



+ 




y.6{x' — x) 



x6{x — x') 



1-x,-kx l-x'.-k'^ 



Figure D.l: The graphs of the effective interaction. Shown are the graphs with an 
effective photon (a), the instantaneous seaguU-graph (b), and the dynamic graph (c) 
and (d) with time ordering (I) and (ff), respectively. 



These functions are listed in Appendix ^ 

Using canonical Hamiltonian matrix elements, as worked out in detail by Tang et 
AL. and listed in Appx. one can easily obtain the same results as using currents. 
As an example, we give the prescription for calculations of effective matrix elements 
originating in the projection of the |ee7)-sector onto the |ee)-sector of the truncated 
Fock-space 

((ee).lKfrl(ee),) = E((ee).|^|(ee7).) , , , , _ , . ((ee7).| V^Kee),), (D.4) 

where V is the vertex operator of a fermion irradiating a photon. For the time ordering 
(I), Fig. (|m])(c), we obtain 



X — x' 
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Appendix D. Calculation of effective matrix elements 



X — x' 



x-x' \ i-x' ^ I -^'i 



The helicity of the irradiated photon is A, in all other respects we use the notation 
of Appx. ly for this vertex operator. The matrix elements for time ordering (II) are 
evaluated in a straightforward way. They evolve from the functions displayed above 
by exchanging particles and anti-particlesQ and a change of the sign because of the 
^-function. We have to multiply these functions by the inverse Hamiltonian, 



((ee7)fc|u; - H\{ee^)k)' 



which is the equivalent of the energy denominator, Eq. ( p.2| ), in our case. Together with 
the much simpler matrix elements of the seagull interaction Fig. (p.l|) (b), we obtain 
the functions of Table ( [bMI ) by substituting all photon helicities by those of fermions. 
This is in accordance with and our previous treatment via currents in Eq. ( p.3|) . 
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^That is, (a;, ki_) ^ (1 - .t, -fc_i). 



Appendix E 
Renormalization 



Because of the divergent graphs in our theory, we have to renormahze the masses of the 
electron and positron, respectively. The singularities are due to the electromagnetic self 
mass of the electron. Fig. ( |2.4| ), which is a diagram of the iterated vertex interaction, 
and the contributions of the two contraction graphs. Fig. (|2.5| )(a+b). The description 
of the renormalization scheme to be presented here follows the line of arguments of 
Krautgartner et al. m, Appx. A]. 



E.l The electron sector 

As a first step, we consider the electron sector, i.e. we restrict the Fock space to contain 
only the states 

|e) = 6i(x=l,fcx=0)|0) 

and 

|e7) = (x, k^)a{^ (1 - x, -A-^) |0) 

which constitute the P- and Q-space analogously to the case of positronium described in 
Chapter 0. We obtain a 2 x 2 block matrix, but here the P-space contains only one state 
for fixed helicity. Using the point symmetries C and H as described in Appx. ^ this 
reduces to effectively one single state in each symmetry sector. Moreover, the Q-sector 
contains only diagonal parts. The seagull graphs are absent due to the properties of the 
theory of effective interactions^. Therefore, the inversion in Q-space becomes trivial, 
yielding 

X l — x' 

The contribution of the iterated interaction is therefore 

^See Chapter || or the plausibiHty argument in Chapter ^. 
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Appendix E. Renormalization 



27r2 Jd ^l-x' 
1 

U L + 

X 1—x' 



mj{l — x') 



i\2 



k' 



X ■ 



{i-x'Y 



1 + -^ 



X ■ 



The domain of integration is restricted by imposing the condition 



Z'2 



X 



+ — — < + ml 

1 rr,' — J 



X' 



This means that the difference of the free invariant masses of the states before and after 
an interaction must not exceed a cutoff A, which plays the role of a regulator. We use 



again the fixing of the parameter oo, Eg. ( 2. 10 ), and arrive at 



W 



a 

2^ 



dx'cfk' 



2m) 



D 



m){l - x'f + k'l x'^ {l-x'Y 



(E.l) 



The evaluation of the contraction terms is achieved in the following way. We extract 
from Table (|C.1|) the contributions of the instantaneous fermion and photon part and 
perform the continuum limit 



C^''\x,k±,X) 
C^^\x,k^,X) 



a 



27r2 Jo 



a 



27r2 



dy 



dy 



oo _, 1 

d^k'^- 

-oo 2 
oo _ 

21:1 



1 



d'k 



yix-y) vix-y) 
1 1 



When we use the principal value as a regulator for these integrals ||5^, we obtain^ 

^ dy I d'^k'j_-, 
J-00 y 

d% 4. 



C^'^\x,k^,X) 



2tt'^ X Jo 
a 



27r2 Jo 



dy 



We have to regulate the integral, since both integrands are singular. If we use here the 
same integration domain as for the iterated interaction, rename the variables and sum 
the two contraction contributions 



C :-- 



/ dx d'^k\ 
27r2 Jd ^ 



{l-x')' 



+ 



1 



X' 



we find that the quadratically divergent term of the instantaneous photon interaction 
is compensated exactly by the analogous expression in the iterated vertex interaction. 



^For the subtleties of these integrals of. |7^, Appx. B]. 
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E.2. The positronium sector 



Eq. Additionally, the fermion part proportional to 1/1 — x' cancels and only a 

logarithmically divergent term survives. The self mass 

Am^ ■.= W + C = / dxd^k'^ I ^ I 

27r2 Jd ^ ym}{2 - x'f + k'l j 

(E.2) 

is therefore only logarithmically divergent, although its three contributions diverge 
quadratically with the cutoff A. 



E.2 The positronium sector 

In the positronium sector as defined in Chapter ||, we have the two states |ee) and {ee^j). 
The calculations of the contractions and the iterated vertex interaction are performed 
exactly as in the previous section. The domain of integration is given by 



m 



j{x - x'y + {xk'_^ - x'k±f 



1 



XX [X — X ) 

The expression for the electromagnetic self mass of the electron (and positron) is 

AM^ = dxd^k'J . - + ^^. 

27r2 Jd ^ \m}{x - x'Y + {xk'^ - x'k^Y xx' x{x - x') ^ 

We can retrieve an expression similar to Eq. ( [E.2|) by substituting 

y :=— and pi := k'^ - yk±, 

X 



where the self masses are related by 

AM^ = -Am^. 

X 

The renormalization is performed by setting 

ruf = rrie = 511 keV, 

or in other words, we adjust the contraction terms so that their divergent pieces cancel 
exactly those of the iterated interaction, i.e. Am^ = 0. 

The renormalization is restricted to P-space. In principle, the counterterms will 
depend on the Fock sector considered. If one takes into account higher Fock states, 
these counterterms will in general not be analytically calculable. They are essentially 
non-p er t urbat ive . 
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Appendix F 

Tables of effective matrix elements 



F.l Introduction and Notations 

The dependence of the effective interaction on the hehcities of in- and out-going particles 
is displayed in the form of tables. A part of the compiled functions has been discussed 



by Krautgartner et al. ||2^. The matrix elements of the effective interaction de- 
pend on the one hand on the momenta of the electron and positron, respectively, and 
on the other hand on their helicities before and after the interaction. The latter depen- 
dencies occur during the calculation of these functions M{x, k±, Ai, A2; x', k'j_, X[, Ag) as 
complicated Kronecker deltas. A good survey is obtained, however, when displaying 
these functions in the form of a table. 

In the remainder of this appendix, the tables are given for the general, angle- 
dependent effective matrix elements. Table ( [l^'.lD , for the matrix elements of arbitrary 
Jz, after integrating over the angles. Table ( [l'\2| ), and for the matrix elements of the 
annihilation graph. Table ( |F.3| ). The following notation is used for functions of the 
type F{x, k±; x', k'jj: 

• An asterisk denotes the permutation of particle and anti-particle. 

F^{x, k±; x', k'j_) := F^i^l — x, —k±; 1 — x', —k'^). 



If the function additionally depends on the component of the total angular mo- 
mentum Jz = n, a tilde symbolizes the operation 

FAn) = F(-n), 
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F.2. General helicity table 



F.2 General helicity table 



final : initial 


(A'i,A'2) =TT 


(A'i,A^) =Ti 


(A'i,A'2) =iT 


(A'i,A^) =ii 


(Ai,A2) =TT 




E*{k,k') 


-E,ik,k') 





(Ai,A2) =Ti 


E*{k',k) 


E,{k,k') 


E^ik.k') 


E,ik',k) 


(Ai,A2) =iT 


-Es{k',k) 


E4{k,k') 


E2{k,k') 


-E;{k',k) 


(Ai,A2) =ii 





Esilk') 


-m,k') 


Ei{k,k') 



Table F.l: General helicity table of the effective interaction. 



A description of the calculation of the functions displayed in the above table was 
given in some detail in Appendix]^. Note that these matrix elements depend in general 
on the vectors k± and k'j_. The functions Ei{k, k') := Ei{x, k±; x', k'j_) read 



Ei{x, k; x', k') 



E4{x, k; x' , k') 



a 1 
2^V 



\xx' {\ — x)\\ — x')) xx'\\ — x){\ — x') 



V \x{\ - x) x'{l-x')J 

Es{x,k;x',k') = -7^,^^(k'^e-^^' -k^\^e-^-^ 

Ztt'^ V xx' \ 1 — X . 



a m'j [x' — xY 
27r2 V xx'{l- x'){l- x)' 



The energy denominator is defined as 



V{x, kx_] x\ k'j_) 



~(x 



J\2 



nil 



+ 



XX' 



- {kl + k'l) + {x- x') 



(1 -x)(l -x' 

47 1 

2 \l-x' 



+ 2k±k'j_ cos(v9 — ip') 



X' 



kl 



1 — X 



X, 
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Appendix F. Tables of effective matrix elements 



F.3 The helicity table for arbitrary 

For an arbitrary = n with n ^ Z one obtains, following the description given in 
Chapter |, the helicity table (Q. 



final : initial 


(A'i,A'2) =Tt 


(a;,a^) =n 


(A'i,A'2) =iT 


(A'i,A'2) =ii 


(Ai,A2) =TT 


Gi(l,2) 


G*(l,2) 


G3(l,2) 





(Ai,A2) =n 


G*(2,l) 


G2(l,2) 


^4(1, 2) 


-^3(2,1) 


(Ai,A2) =iT 


^3(2, 1) 


G4(l,2) 


G2(l,2) 


-Gl{2, 1) 


(Ai,A2) =ii 





-^3(1,2) 


-G|(l,2) 


Gi(l,2) 



Table F.2: Helicity table of the effective interaction for = ±n, x > x' . 



Here, the functions ^,(1,2) := Gj( fc^) are given by 

'1 1 ~ 



^ ^ xx' ' (1 — x)(l — x') 
kik' 



Int{\l-n\) 



G2(x, k±; x', k\ 



xx'il — — x') 

nil 1 



-Int{\n\) 



+ 



+ ^^^+ Int{\n\) 



^ \xx' ' {1 — x){l — x') ) ' x{l — x) ' x'{l — x') 
+k^k\ 



Int{\l - n\) _^ Int{\l + n\) 



XX' 



(1 -x)(l -x') 



G'3(x, x', k'^ 



-ruf — - 



1 — x' 

k'Jnt{\l - n|) - k± Int{\n\) 

I — X 



Gii^x, k±] x' , k'j_) = — 



(X — X 



/\2 



Int{\n\). 



^ xx'{l — x'){l — x) 
The function Int{n) is defined as 

Involved in this expression are the definitions 



Tfl I X 

[X-x')'^—f I : + 



1 



--[x-x 



2 \xx' 



1 — x' x' 



- kl + k'l 

1 1 



1 — X X 



92 



F.4. Helicity table of the annihilation graph 



and 

A ' 

/a2 - 4A;2 



B = i(l-aA) 



F.4 Helicity table of the annihilation graph 



finahinitial 


(A'i,A'2) =TT 


(A'i,A'2) =n 


(Ai,A^) =iT 


(A'i,A'2) =ii 


(Ai,A2) =TT 


Fiil,2) 


i^3(2, 1) 


F3*(2,l) 





(Ai,A2) =n 




f;{i,2) 


i^4(2, 1) 





(Ai,A2) =iT 




F4(l,2) 


i^2(l,2) 





(Ai,A2) =ii 















Table F.3: Helicity table of the annihilation graph for > 0. 



The functions 2) := Fi{x,k^;x',k'j_) were calculated in Chapter 



Fi{x,k±;x',k'j_) 
F2{x,k±;x',k'j_) 
Fs{x, k±;x', k'i_) 
Fa{x, k±;x', k'j_) 



a 2rn? ( 1 



X 



xj\x' 1 — x' 



|Jb|,1 



a 

TT 



2 k^k'j 



U XX 



k\ 



a 2m / 1 1 
^Ai - + - 

TT iO* \X \ — X J \ — X 



a 

TT 



2 kik\ 



UJ* x'{l — x) 



J. 1,1 



45 J., 



The table for = — 1 is obtained by inverting all helicities. Note that the table has 
non- vanishing matrix elements for | J^l < 1 only. This restriction is due to the angular 
momentum of the photon. 
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Appendix G 
Numerics 



To solve the integral equation ( |2.11| ), central to this work, we apply several numerical 
methods, which are described in the sequel. 

G.l Transformation of the variables 

The invariant front form wavefunctions coincide with the equal-time wavefunctions 
defined in the infinite momentum frame p. 379]. Let the relative momentum of the 
i-th particle in a composite system be Pi. In the center of mass system holds 

Y.P^ = 0. (G.l) 

i 

We want to express the variables x, in terms of the equal-time variables p±,Pz- With 



(|G.1| ) we have 



pt = E + pz 



vt = E -pz 



where E = ^Jvrfi + p^. If we write the relative momentum in polar coordinates 

p={fi sin 9 cos ip, fismO simp, fi cos 9) , 
it is clear from the definition of the front form momentum fractions ( |A.1| ) that 



X 



1 / LL cos 6 

2'^ + 



k± = li± = {fisinO cosip, ^sm9 simp). 
The inverse transformation is 



cos 9 



k'i + m?{2x- ly 
\ 1 - (2x - 1)2 

-2x) 



\ kl + m^{2x-iy' 



(G.2) 
(G.3) 
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G.2. Symmetries 



These are exactly the coordinates as used by Karmanov for the deuteron wave- 
function and by Sawicki [|6|, Eqs. (2.3), (2.4)]. The Jacobian reads 

, 1 + u2(l - cos2 6') 9 
J = ] 3 -n''sm9. 

For the integration measure holds 



dx d^k± = I dfi dcos6 / d(f ■ 

Jo J-OD Jo J-1 Jo 



2 (m2 + /i2)3/2 • 
The new variable n is referred to as as off-shell mass, because of 

The restriction on the momenta ( p.l3| ) translates into 

4(mJ + /i^) < a2 + (G.4) 

and is therefore a restriction of the off-shell mass. 

In the next paragraph the integration over the angular variable ip is described. After 
this, only two variables are left in the game: /z and cos^. The integrations over these 
coordinates are discretized, i.e. the integral is mapped onto a sum, with help of the 
Gauss-Legendre algorithm (cf. e.g. |^). The domains of the integrations are chosen 
in such a way, that 

/^e[o,-], cos0e[-i,i]. 

Here, abscissae and corresponding weights are denoted by 

fii,Wi i = l..Ni 
cos 9j, Wj j = I..N2. 

With A — i> cx) we have to map the interval [0, 00] for numerical integration onto a finite 
domain. We follow [|4|, p. 60, case 7] and use as a mapping function 

1 + /i 

The mapping of the weights is calculated straightforwardly. 



G.2 Symmetries 



First, we want to restore the symmetry of the Hamiltonian matrix broken by the above 
transformation of the variables. In discrete variables, the transformation has the Jaco- 
bian 

^ m^ + /xf(l-cosgj) 

ttij :- fi^ —j=== . [^.t)) 

mj + fif 
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Appendix G. Numerics 



From this transformation an unsymmetric matrix evolves. As a consequence, the nu- 
merical effort to diagonalize it increases considerably. We symmetrize the matrix there- 
fore by substituting for the wavefunction 

The so-modified matrix elements of the effective interaction are 



{x, k±; Ai, \2\Ves\x', k'^; \[, \'2)^JwiWjW'^w'iaija'f.^. 

Now, we address to the implementation of the point symmetries of the theory into 
the matrix elements. The Lagrangian density of QED is invariant under the trans- 
formations of time reversal T, parity V, and charge conjugation C. Opposed to that, 
the Hamiltonian in front form dynamics is invariant under the latter operation only. 
The other operations can be used to construct further, in general more complicated, 
symmetries. 

We consider as a further symmetry of our Hamiltonian matrix the charge conju- 
gation, which is represented mathematically by the unitary operator Uc- Eigenstates 
to this transformation have the eigenvalues ttc = ±1 and the creation operators of the 
different particles transform as follows ||9^ 

Ucb[{x,k^)Uc' = d[{x,k^), 
Ucd{ix,k^)Uc' = h{{xM), 
Uca{{x,k_i_Uc^ = -a{{x,k±). 

The helicities are A = ±1. We construct out of arbitrary P-space wavefunctions with 
undefined charge 

l^P) = -^eeix^k^; Xi, X2)b\^ix,k^)d\^{l - X, -k^_)\0) 
eigenstates to the charge conjugation operator 

kc = ±1) = 

^ (G.6) 

-j=^ee{x,kr,\i,\2) [b{^{x,k±)d{^{l - X, -k±) T bi^i^ - X, -k_L)d{^{x,k_L^ |0). 

The remark is in order, that if the quantum numbers of the creation operators coincide, 
there is only one eigenstate, having the eigenvalue ttc = —1. This explains the different 
dimensions of the blocks of the Hamiltonian matrix, which are build of eigenstates to 
the point symmetries. 

As mentioned earlier, further invariants can be constructed from the two remaining 
symmetries of the Lagrangian density [0. We use 
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G.2. Symmetries 



In Chapter ^ of the present work we address to the special case of Jz = 0. In this 
case, we can use the combination Ti. := VT as another symmetry. We give again 
the behaviour of the creation operators under the transformation via the anti-unitary 
operator[] Vn 

Physically spoken, this operation inverts the helicities and causes a phase change, de- 
pendent on the choice of the spinors ||9^[0- One can read off the eigenfunctions of 
the operator Vn 

|7r^ = ±l) = ■^[^ee{x,kr,\i,^2)b{^{x,k^)d{^{l-x,-k±) (G.7) 

±(-l)^^+^^^^,l^{x, kr, Ai, A2)&U k^)d-xA^ - -^^)] |0)- 

Applying the two point symmetries, the whole Hilbert space is decomposed into 
four subspaces with well-defined charge and "T-parity", and the Hamiltonian matrix 
becomes block-diagonal. As a consequence, only four much smaller matrices need to 
be diagonalized. The numerical effort decreases with the third power of the matrix 
dimensions, and a lot of computer time can be saved. The construction of the simul- 
taneous eigenstates to C and 7i is supplied by the computer routine GENETIX, which 
generates the whole Hamiltonian matrix by using a subroutine, which calculates the 
unsymmetrized matrix elements. 

Finally, we use for the implementation of one further symmetry the fact, that the 
Lagrangian density and the Hamiltonian operator of QED are invariant under rotations 
in the x-y-plane in front form dynamics, too. Due to this, we can introduce, as another 
quantum number, the 2;-component of the total angular momentum J^, to classify our 
states. The angles enter via 

_ f k± cos ip \ 
\ k±sm(p J ' 

The new effective interaction is obtained by Fourier-transforming the old one 

{x, k±, Lz] Ai, X2\Ves\x\ k'j_, L'/, X[, X'^) : = 

^ r e-^^^'^ r d^' e^^^^'(x,fc^,<^;Ai,A2|Kfr|2;',fc^,<^';A;,A^). 

But, neither nor Sz = Xi + A2 are good quantum numbers and we set Lz — Jz — 
5*^. From these calculations we gain the functions of the helicity table (|F.2|) . This 
functions are used in the computer code in the subroutine PHYSIX, which calculates the 
unsymmetrized matrix elements. 

^It is the direct product of the unitary operator of the parity transformation and of the anti-unitary 
operator of the time reversal transformation Vn := Vp (B Vr ■ 
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Appendix G. Numerics 



G.3 Coulomb trick 

The underlying structure faced in treating bound states in Quantum Electrodynamics 
is the Coulomb problem. Obviously, in the non-relativistic limit, the derived effective 
integral equation ( |2.11| ) will be mapped into the Coulomb-Schrddinger equation 0. 

In the introduction, we have already made the point that momentum coordinates are 
essential, because they render the exact treatment of the center of mass motion possible. 
A major part of the numerical problems in solving Eq. (p.ll|) is corresponding to the 
solution of Coulomb problems in momentum space. This problem was examined in 



detail in |^4[ . We present the numerical techniques applied to the treatment of integral 
equations with a weakly singular kernel. 

In the Coulomb Schrodinger equation a kernel is involved with a point singularity 
which is analytically integrable. Nevertheless, the matrix elements become senseless at 



that very point. Therefore, we apply the so-called CoULOMB trick [^, which is known 
as the "Nystr0m method" in the mathematical literature. 

In momentum representation the Schrodinger equation of the hydrogen atom reads 

- ^ fd'p'-^ = Em. 

2m^^^' 2n^J ^(p'-p')2 ^^^^ 

For simplicity, we consider only S-waves, i.e. rotationally invariant states. We can 
integrate out the angular variables and using Bohr units we arrive at 

+ -(dp'^- In ^(P') = ^^(P)- 

-K J p [[P + Pj ) 

Next, we translate the integral with help of the Gauss-Legendre algorithm in a sum 



r- 



" Pi \{pi+Pjy 



This discretized equation has a singularity at Pi = pj. To treat it, we add and subtract 
two terms which cancel exactly in the continuum limit. Both are diagonal, since the 
divergence occurs on the diagonal only. The new equation looks like 

P.Vfe) + - E^.^ 1^ ir^T^l t^^^^) - 9iP^,PMP^)] (G.8) 
^ j=l Pi I \Pi ' Pj ) ) 

1 f J lP\ \iPi~ P'^'^ 



+- f dp'P- In W) = Ei^ipd- 

TT Pi \ [Pi + p'y J 



One convinces oneself easily, that the second term in the squared bracket, multiplied 
by a function g{pi,Pj) generating better convergence, nullifies this bracket as pi = pj. 



^For a proof see [[73| or [gg Appendix C]. 
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G.3. Coulomb trick 



For oo it is compensated exactly by the additional integral expression. It must 

hold 

g{p,p) = 1. 

Since the ground state wavefunction is known, one sets for instance 

The integral in ( |G.8D , the continuous part of the Coulomb trick, will in general not be 
solvable analytically. It suffices instead, to evaluate it numerically with a much higher 
precision than achieved by doing the sum, i.e. the discrete part. Here, we obtain the 
function 

{(p. — p I 
^p.+p.)2 j [^fe) - 9{P^,Pj)■ip{P^)] ■ 

It can be continued continuously into p = p' and therefore is integrable numerically 
much easier than the original function. There are even possibilities of improving this 
method |9^. The matrix is no longer symmetric. We substitute 

(piPi) '■= VWiPi'ipiPi) 

and obtain finally 

pf4>(pi) H — JWiWj In < ^— \ Mpj) + ( Coulomb- Counterterms) = E(j)ipi). 

T^J^l {{Pi+Pj) ] 

The result of the manipulations is a considerable improvement of convergence, docu- 
mented in Ref. 
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Appendix H 

Listing of the computer code 



To render the verification of tlie results of the present work possible and to demonstrate 
the implementation of the techniques described throughout the text in praxi, we show 
the listing of the program, written in C. 



H.l Description of the program 

The program MesonixQ is constructed in a modular fashion, i.e. one problem corre- 
sponds to one subroutine. The whole program consists of the below parts or subrou- 
tines: 



meta_mesonix: This part of the program was written to render its application 
more convenient for the user. An input file for the main program mesonix is 
created. It contains all parameters, such as cutoff, fermion mass, number of 
integration points, etc. Furthermore, small program tools are supplied, which 
enable subsequent calls of mesonix to be handled directly. The function vary_J_z, 
for example, was created to calculate the spectra for different values of Jz at a 
fixed discretization. 

mesonix: The main program calls the subroutines according to the parameters 
of the input file input .dat. 

numerix: Using the NAGLiB- Routine dOlbcf [Q, the abscissae and weights for 



the numerical integrations are obtained using the Gauss-Legendre algorithm. 

coulomb_trix: Here, the calculation and initialization of the counterterms of 
the Coulomb trick are processed. It decides whether the results of a former 
calculation are read in from a file, or if they are (re-) calculated for a different set 
of parameters. 



"'^ Mnemonic for mesonzcs= "anything, concerning mesons" , like numerics= "anything, concerning the 
numeric aspects" . Here, a meson is defined as a fermion-antifermion system. 
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genet ix: This routine governs the construction of the symmetrized Hamiltonian 
matrix. By using the charge conjugation symmetry C and the combined parity and 
time reversal symmetry Ti, four matrices which describe the interaction between 
the eigenstates of these symmetries, are generated. 

asymmetrix Analogously, this procedure governs the construction of the unsym- 
metrized Hamiltonian matrix. In the case J^t^O the parity and time reversal 
operation is not a symmetry any more. A special treatment is necessary, which 
increases the numerical effort enormously. 

physix: The matrix elements are calculated by evaluating the function listed 
in the helicity tables of Appendix 0. The Coulomb counterterms are introduced 
when all quantum numbers of the incoming and the outgoing states are identical. 



• arithmetix: With help of the NAGLiB routine f02abf |^9|, the Hamiltonian 
matrix is diagonalized and the spectrum, as well as the wavefunctions are calcu- 
lated and sorted (mOlcaf). It is possible to store the wavefunctions. For this, 
they are translated into the original normalization. 

• publix: This routine formats and outputs the results. 

All global variables are explained in the listing. The program is written in standard 
C, but makes no excessive use of typical C structures and should consequently be 
understandable to readers familiar with Fortran. 
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Meta_mesonix 



/* */ 

/* Program : META.MESONIX */ 

/* */ 

/* Class : Main program */ 

/* */ 

/* Purpose : META.MESONIX calls MESONIX and manipulates the */ 

/* input file to vary different outer parameters */ 

/* via INPUT_MES.DAT. */ 

/* */ 



/************************************************************************/ 



#include <stdio.h> 
#include <math . h> 



double alpha; 

char default .directory [80] ="OUTPUT/" ; 
char special_id[80]="_"; 
char *name; 
FILE *fp; 



void create_file (double beta, double lambda, int nl, int n2, int number, 

int J_z, int trix) 
/* create input file for 'mesonix.c' */ 
{ 

double m; 

int test,Anni,print_EF; 

char *name_in = "Input_mes.dat"; 

FILE *fp; 

m =1.0; /* fermion mass */ 

test =1; /* test variable */ 

Anni =0; /* control for incl . /excl . of annihilation graph */ 

print_EF =0; /* control for writing the eigenf unctions in file */ 



if (J_z!=0) test += 4096; /* asymmetric case */ 



fp = fopen(: 
f printf (f p, 
fprintf (fp, 
f printf (fp, 
fprintf (fp, 
fprintf (fp, 
fprintf (fp, 
fprintf (fp, 
fprintf (fp, 
fprintf (fp, 
fprintf (fp, 
fprintf (fp, 
fprintf (f p, 
fprintf (f p, 
fprintf (f p, 
fprintf (f p, 
fprintf (f p, 
fprintf (f p, 
fprintf (f p, 
fprintf (f p, 
fprintf (f p, 
fprintf (fp, 
fprintf (fp, 
fprintf (fp, 
fprintf (fp, 
fprintf (fp. 



naine_in, "w") ; 



MESONIX.input.f ile\n") ; 



-\n"); 
-\n\n"); 



f ermion_mass=\n") ; 

"/.18.12f\n",m) ; 
cut-off=\n") ; 

"/.18.12f\n", lambda) ; 
coupling=\n") ; 

"/.18.12f\n",beta); 
test=\n") ; 

'/.3d\n" .test) ; 
J_z=\n") ; 

■/.3d\n",J_z) ; 
Anni=\n") ; 

'/.3d\n" ,Anni) ; 
Nl=\n") ; 

■/.3d\n",nl); 
N2=\n") ; 

■/.3d\n",n2); 
number _EV=\n") ; 

■/,3d\n" .number) ; 
trix=\n"); 

•/.3d\n",trix); 
def ault_directory=\n") ; 

%s\n" .default .directory) ; 
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fprintf (fp, "special_id=\n") ; 
f prlntf (f p, " 7.s\n" , special_id) ; 
fprintf (fp, "print_EF=\n") ; 
fprintf (fp, " "/.SdNn" ,print_EF) ; 

fprintf (fp," \n"); 

f close (f p) ; 

} 

void vary_coupling(int nn,int N) 

/* vary the coupling ALPHA: Nl=N2=nn, */ 

{ 

int 1 , j , k ; 
double alpha ; 
FILE *f_in,*f_out; 

char *name_out = "EVs_couplinx.dat"; 

f_out = fopen(name_out, "w") ; 

f close (f_out) ; 

for (1=1; i<=N; ++i) 

{ 

f_in = fopen(name,"w") ; /* create output file for MESOHIX */ 
f close (f_in) ; 
alpha = 10.0*i/N; 

printf ("META_MESONIX: ALPHA = 7.8. 5f\n\n" .alpha) ; 

Great e_f lie (alpha, 1.0,nn,nn,10,0,l); 

systemC'mesonix.out") ; 

f_in = f open (name, "r") ; 

f_out = f open(name_out , "a") ; 

f scanf (f _in, '"/.d" ,&k) ; 

fprintf ("k= y.3d\n",k); 

fprintf (f.out, "'/.IS. 12f\n" .alpha) ; 

for (j=l; j<=10; ++j) 

{ 

f scanf (f _in, '"/.If" ,&alpha) ; 
printf ("EW = '/.IS . 12f\n" .alpha) ; 
fprintf (f .out , ""/.18 . 12f \n" . alpha) ; 

} 

f close (f_in) ; 
f close (f_out) ; 

} 

} 

void vary_N(int begin_step, int end_step, int J_z) 

/* vary number of mesh points from <begin_step> to <end_step> */ 

{ 

int i ; 

sprintf (name , "'/,sEigenValues/EVs'/,sJ'/,d. dat " . 

def ault_directory .special_id, J_z) ; 
printf ("y,s\n" .name) ; 

fp = fopen(name."w") ; /* create output file for MESONIX */ 

f close (fp) ; 

for (i=begin_step; i<=end_step; i+=2) 
{ 

printf ("\n\nMETA_MESONIX: Nl = '/.Sd. N2 = "/.SdNn" . i . i) ; 
create.f lie (alpha. 1 . . i , i . 500 , J_z . 0) ; 
system ( "mesonix . out " ) ; 

} 

} 

void special_case (double lambda, int nl, int n2, int J_z) 

/* different Nl. N2 */ 

{ 
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int i ; 

spr intf (name , "'/,sEigenValues/EVs'/.s J%d . dat " , 
default .directory , special_id , J_z) ; 
printf ("'/,s\n" ,name) ; 

fp = f open(name, "w") ; /* create output file for MESONIX */ 

f close (f p) ; 

printf ("\n\nMETA_MESONIX: Nl = 7,3d, N2 = '/.SdXn" ,nl ,n2) ; 
create_f ile (alpha , lambda , nl , n2 , 500 , J_z , 0) ; 
systemC'mesonix.out") ; 

} 

void Coulomb_J_z(int end_step,int max_n) 
/* calculate Coulomb counterterms */ 
{ 

int i ; 

for (i=0; i<=max_n; ++i) 
{ 

printf ("\n\nMETA_MESONIX: J_z = •/.3d\n",i); 
sprintf (name , "'/,sEigenValues/EVs7,sJ%d.dat" , 

def ault_directory , special_id , i) ; 
fp = f open(name, "w") ; /* create output file for MESONIX */ 

f close (f p) ; 

create_f ile(alplia, 1 . ,end_step, end_step,500, i , 1) ; 
system ( "mesonix . out " ) ; 
if (i>l) 
{ 

printf ("\nMETA_MESONIX: J_z = •/.3d\n",i); 
sprintf (name , "%sEigenValues/EVs'/,sJ'/,d. dat" , 
def ault_directory ,special_id,-i) ; 

fp = f open (name, "w") ; /* create output file for MESOKIX */ 

f close (fp) ; 

create.f ile (alpha, 1.0, end.step , end.step , 500 , -i , 1) ; 
system ( "mesonix . out " ) ; 

} 

} 

} 

void special_Coulomb_J_z(int Nl, int N2, int J_z) 

/* calculate Coulomb counterterms for special case: (Nl,N2,J_z) */ 
{ 

printf ("\n\nMETA_MESQNIX: J_z = "/.SdXn" , J_z) ; 
create_f ile (alpha , 1 . , Nl , N2 , 500 , J_z , 1) ; 
system("mesonix.out") ; 

} 

void vary_J_z(int begin_N, int end_N, int max_n) 

/* vary J_z from -max_n till +max_n */ 

{ 

int i ; 

for (i=0; i<=max_n; ++i) 
{ 

printf ("\n\nMETA_MESONIX: J_z = •/.3d\n",i); 
vary_N (begin_N , end_N , i ) ; 
if (i>0) 
{ 

printf ("META_MESONIX: J_z = '/.SdXnXn" ,-i) ; 
vary _N (b6gin_N , end_N , - i ) ; 

} 

} 
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} 



double LambdaCint scale, Int n) 

/* returns Lambda in different scales , n=l . . 10 */ 
{ 

switch (scale) 
{ 

case 0: return n*l . 0*6 . 0*alpha; break; /* Krautgaertner scale */ 
case 1: return n; break; /* linear scale */ 

case 2: return (15 . 0+5 . 0*n) ; /* linear till 65 scale*/ 

break; 

} 

return 1.0; /* default */ 

} 



void create_environement() 

/* create directories for output */ 

{ 

FILE *fp; 



printf ("\nMETA-MESONIX \n\n") ; 

fp = fopen("create_dir","w") ; /* create file testing if dir. exists */ 
fprintf (fp,"if [ ! -d $1 ] \nthen mkdir $l\n"); 
fprlntf (fp, "echo \"create_dir : Directory '$1' created. \"\n") ; 
fprintf (fp, "else echo \"create_dir : Directory '$1' exists!\" \nfi"); 
f close (f p) ; 

naine=(char *)malloc(80) ; 

sprlntf (name , "/bin/bash create_dir °/,s" ,def ault_directory) ; 
system (name) ; 

sprintf (name , "/bin/bash create_dir 7,s7,s" ,def ault_directory , 

"Eigenvalues/") ; 

system (name) ; 

sprintf (name , "/bin/bash create_dir '/.s'/.s" ,def ault_directory , 

"EigenFunctions/") ; 

system (name) ; 

sprintf (name , "/bin/bash create_dir °/,s7,s" ,def ault .directory , 

"Coulomb.Trick/") ; 

system (name) ; 

system("rm -f create_dir") ; /* remove batch file */ 

sprintf (name, "/ismesonix.log" ,def ault_directory) ; 
printf ("Creating ''/,s' . . .\n",name) ; 

fp = fopen(name,"M") ; /* create LOG file for MESOKIX */ 

f close (fp) ; 

printf (" \n\n") ; 



void vary_Lambda(int J_z) 
/* vary cutoff Lambda */ 
{ 

int i,nl,n2; 
double cutoff; 



ni = 41; 
n2 = 11; 

for (i=l; i<=5; ++i) 
{ 

cutoff=Lambda(0,i) ; 

printf ("\n\nMETA_MESONIX: Lambda = '/.lO .8f\n" , cutoff ) ; 
sprintf (special_id, "_L'/,ld_" ,i-l) ; 
printf ("7,s\n" ,special_id) ; 

create_file(alpha, cutoff ,nl ,n2, 500, J_z, 1) ; /* Coulomb terms */ 

system("mesonix.out") ; 

special_case(cutoff ,nl,n2,J_z) ; /* N1!=N2 */ 

} 

for (i=6; i<=10; ++i) 
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{ 

cutoff=Lainbda(0,i) ; 

printf ("\n\nMETA_MESONIX: Lambda = 7,10.8f\n" .cutoff ) ; 
sprintf (special_id, "_L'/,ld_" , i-1) ; 
printf ("'/,s\n" , special_id) ; 

create_file(alpha, cutoff ,nl ,n2, 500, J_z, 1) ; /* Coulomb terms */ 
systemC'mesonlx.out") ; 

special_case(cutoff ,nl,n2, J_z) ; /* N1!=N2 */ 

} 

for (i=l; i<=10; 
{ 

cutoff =Lambda(2 , i) ; 

printf ("\n\nMETA_MESONIX: Lambda = 7.10.8f\n" .cutoff ) ; 
sprintf (special_id, "_LXy,ld_" , 1-1) ; 
printf ("'/,s\n" , special_id) ; 

create_file(alpha. cutoff .nl.n2, 500, J_z,l) ; /* Coulomb terms */ 
system ( "mesonix . out " ) ; 

special_case(cutoff ,nl,n2, J_z) ; /* N1!=N2 */ 

} 



void Lambda20(int J_z) 

/* calculate for Lambda=20 m_f . Nl=41, N2=ll */ 
{ 

int nl=41; 
int n2=ll; 
double cutoff =20.0; 



sprintf (special_id, "_120n4711_") ; 
printf ("7,s\n" , speclal_id) ; 

create_file(alpha, cutoff ,nl ,n2. 500, J_z. 1) ; /* Coulomb terms */ 
systemC "mesonix . out") ; 

special_case(cutoff .nl.n2. J_z) ; /* N1!=N2 */ 



void mainO 
{ 

alpha = 1.0/137.0359895; 
alpha = 0.3; 



/* coupling constant */ 



create_environement () 
Coulomb. J_z (21,3); 
vary_N(5,21,0) 
vary_N(5,21,l) 
vary_N(5,21,2) 
vary_N(5,21,3) 

vary_Lambda(0) ; 

} 



/* create directories for output */ 

/* calc. counterterms */ 

/* calc spectra as func. of H. J_z=0 



/* calculate spec, for Lambda= l...,55 m_f*/ 



.,3*/ 
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Mesonix 



/* */ 

/* Program :MESONIX */ 

/* */ 

/* Class : Main program */ 

/* */ 

/* Purpose : MESONIX calculates the eigenstates and -values */ 

/* of a system consisting of two particles with */ 

/* same mass, i.e. of mesons. */ 

/* */ 
/************************************************************************/ 

/************************************************************************/ 

/* */ 

/* Structure of the code : MESONIX */ 

/* =============== */ 

/* */ 

/* MAIN PRG SUBROUTINES SUB-SUBROUTIKES */ 

/* */ 

/* */ 

/* main > NUMERIX */ 

/* I */ 

/* I > COULOMB.TRIX */ 

/* I */ 

/* I > GENETIX > PHYSIX */ 

/* I */ 

/* I > ASYMMETRIX > PHYSIX */ 

/* I */ 

/* I > CALCULATIX */ 

/* I */ 

/* I > PUBLIX */ 

/* */ 

/* */ 

/* */ 

/* Description of the subroutines : */ 

/* */ 

/* */ 

/* NUMERIX calculates the weights and abscissae for the */ 

/* numerical integration */ 

/* */ 

/* COULOMB_TRIX calculates/initializes the Coulomb counter terms */ 

/* */ 

/* GENETIX generates the Hamilton matrix */ 

/* */ 

/* PHYSIX provides the single matrix elements */ 

/* */ 

/* CALCULATIX diagonalizes the matrix and sorts eigenvalues */ 

/* */ 

/* PUBLIX controls the output of data */ 

/* */ 

/* */ 

/* */ 

/* */ 

/* Outer Parameters: */ 

/* */ 

/* */ 

/* J_z = n ; z component of total angular momentum */ 

/* */ 

/* Anni = 0/1 ; omit/include annihilation graph */ 

/* */ 

/* trix = 0/1 ; read in from file / calculate Coulomb term */ 

/* */ 

/* number.EV = m ; write m eigenvalues in file 'EVs_Jn.dat' */ 
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print _EF 



0/1 



do not write/do write eigenfimctions in file 
'EFunc_X.dat' 



/* 
/* 
/* 
/* 

/* 

/* 

/* 

/* Test Facilities of 


/* 

/* bit condition 

/* 

I* 

I* all test 
/* 

/* (test & 1) != 
/* 
/* 
/* 

/* 1 (test & 2)!= 
/* 
/* 

/* 2 (test k 4) ! = 
/* 

/* 

I* 

I*Y1 (test &4096)!=0 Supress the symmetries, use ASYMMETRIX 
/* 



*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 
*/ 



the Program : 



action 



== nothing will be tested 



NUMERIX will write the abscissae and weights of 
the numerical integration into the file 
' gauss.vals . dat ' 

GENETIX will write a off-diagonal spin table 
in "SPIN-TABLE_MES.DAT" 

PUBLIX will write the eigenvalues of the block 
matrices in EV_1_MES.DAT, EV_2_MES.DAT, 
EV_3_MES.DAT, EV_1_MES.DAT respectively. 



#include <malloc.h> 
#include <stdlo.h> 
#include <math.h> 
#include <unistd.h> 
#include <signal.h> 

#define PI 3.141592653589793 

#define NNl 41 

#define NN2 21 

#define np 1 

#define Ip 4000 

#define lip 2000 

#define N.sym.l 2*(N1*N2)+2*N1 

#define N_sym_2 2*(N1*N2)+2*N1 

#define nagtest 1 

tdefine old.CT 



/* max. # of integr .points (mu) */ 
/* max. # of integr .points (theta) */ 
/* constant for 'dOlalf, etc. */ 
/* constant for 'dOlalf, etc. */ 
/* constant for 'dOlalf, etc. */ 
/* dimensions of blocks (syrmnetrix)*/ 

/* use NAGLib for diag if !=0 */ 
/* use old(=Krautgaertner) Coulomb */ 
/* trick if != */ 



/* ===================> external functions <=================== */ 

extern long timeO; 

/* time measurements */ 

extern double f02abf_( double *A, int *IA, int *NP, double *R, double *V, 

int *IV, double *E, int *IFAIL) ; 
/* diagonalization routine from NAGLib*/ 

extern double Coulomb_trick_function(double *) ; 

/* Continous Coulomb counterterm, dependend on (mu, theta) */ 

extern double Coulomb_trick_integrand (double *) ; 

/* Continous Coulomb counterterm, integrated over (theta), dependent on (mu) */ 
extern double Coulomb_discrete(int jm, int si, int s2) ; 
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/* Discrete Coulomb counterterm, dependent! on spins and mu[jm] */ 
extern double integrand (double *) ; 

/* Krautgartner continous Coulomb counterterm, analytically integrated over */ 
/* (theta) , dependent on (mu) */ 

extern double old_continous(int jm, int jt) ; 

/* Krautgartner continous Coulomb counterterm */ 



/* 



global variables 



*/ 



double m, lambda, ALPHA, p_bohr 
double aa, A,B, V_seag; 
double x[3] ,k[3] ; 
double *mu,*wmu; 
double *theta,*wtheta; 
double *ma,*r,*v,*e; 
double *mat_l,*r_l,*v_l,*e_l 
double *mat_2,*r_2,*v_2,*e_2 
double *mat_3,*r_3,*v_3,*e_3 
double *mat_4,*r_4,*v_4,*e_4 
double mu_global; 
double CT [4] [NNl] [NN2] ; 

int J_z, test ,N1 ,N2 , steps ,Anni ; 
int trix,number_EV, print _EF; 
int EV_number ; 
int IndexMatrix; 
int IndexMat [4] ; 
int jjl,jj2; 

int spins_parallel,spin_l; 

int J_z_Coulomb ; 

long timel,time2; 

char *order; 

char *name_EV; 

char *name_log; 

char *default .directory [80] ; 

char *special_id[80] ; 

char *CT_name [4] ; 

FILE *f_log; 



/* mass, cut-off, coupling, Bohr momentum */ 

/* variables used in 'physix' and 'genetix' */ 

/* old variables (x,k) instead (mu, theta) */ 

/* mu: abscissae and weights for Gauss-Legendre */ 

/* cos (theta): " " Gauss-Legendre */ 

/* variables for NAGLIB diagonalization */ 

/* variables for MGLIB diagonalization */ 

/* variables for NAGLIB diagonalization */ 

/* variables for NAGLIB diagonalization */ 

/* variables for NAGLIB diagonalization */ 

/* Coulomb trick: global variable for C_T_f unction*/ 

/* matrix for Coulomb data, first index is */ 

/* n = uu,ud,du,dd */ 

/* OUTER parameters , read in from file */ 
/* OUTER parameters , read in from file */ 
/* arithmetix: number of eigenf unction */ 
/* arithmetix. c */ 
/* arithmetix. c */ 

/* global variables for Coulomb trick */ 

/* Coulomb trick: choose function Gl <-> G2 */ 

/* auxiliary variable for Coulomb trick */ 

/* var. for time measurements */ 

/* for NAGLib routine mOlcaf */ 

/* name of file containing the eigenvalues */ 

/* name of LOG file */ 

/* read in by logix, def . dir. for output */ 
/* read in by logix, special ident. for output*/ 
/* names for the Coulomb trick files */ 
/* File variable for 'log' file */ 



/* 



=> functions for error handling 



*/ 



void nrerror(char error_text [] ) 
/* standard error handler */ 
{ 

fprlntf (stderr, "Numerical Recipes run-time error .. .\n") ; 
f prlntf (stderr , "'/,s\n" , error_text) ; 
fprintf (stderr, ".. .now exiting to system. . .\n") ; 
exit(l) ; 

} 



/* 



=> functions for (de) allocation 



*/ 



double *dvector(long nl, long nh) 

/* allocate a double vector with subscript range v[nl..nh] */ 
{ 

double *v; 



v= (double *)malloc((size_t) ((nh-nl+l+l)*sizeof (double))) ; 
if (!v) nrerror ("allocation failure in dvectorO"); 
return v-nl+1: 
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void free.dvector (double *v, long nl, long nh) 

/* free a double vector allocated with dvectorO */ 

{ 

free((char*) (v+nl-1)); 

} 

/* ===================> inclusion of subroutines <====================== */ 

#include "numerix.c" 
#include "physix.c" 
#include "genetix.c" 
#include "coulomb_trix.c" 
#include "arithmetix.c" 
#include "publix.c" 
#include "asynmetrix.c" 

/* ==============> functions used for initalization <=================== */ 



void init_ma(int dim) 

/* initialize matrix for asymmetic NAGLib diagonalization */ 
{ 

ma = dvector (0,dim*dim - 1); 

e = dvectorCO, dim-1) ; 

r = dvectorCO, dim-1); 

V = dvectorCO, dim*dim - 1); 

IndexMatrix = 0; 

} 

void init_matCint diml, int dim2, int dimS, int diin4) 

/* initialize matrices 1-4 for symmetic NAGLib diagonalization */ 

{ 

int i ; 

for Ci=0; i<=3; ++i) IndexMat[i] = 0; 

mat_l = dvector CO, diml*diml - 1); 

e_l = dvectorCO, diml-1) ; 

r_l = dvectorCO, diml-1); 

v_l = dvectorCO, diml*diml - 1); 

mat_2 = dvector CO, dim2*dim2 - 1); 

e_2 = dvectorCO, dim2-l) ; 

r_2 = dvectorCO, dim2-l) ; 

v_2 = dvectorCO, dim2*dim2 - 1); 

if Cdim3>0) 
{ 



mat_3 


= dvectorCO, dim3*dim3 - 


1); 


e_3 


= dvectorCO, dim3-l) ; 




r_3 


= dvectorCO, dim3-l) ; 




v_3 


= dvectorCO, dim3*dim3 


- 1); 


mat_4 


= dvector CO, dim4*dim4 - 


1); 


e_4 


= dvectorCO, dim4-l) ; 




r_4 


= dvectorCO, dim4-l) ; 




v_4 


= dvectorCO, dim4*dim4 


- 1); 



} 

} 

void free_matCint diml, int dim2, int dim3, int dim4) 

/* deallocate matrices 1-4 used for symmetric NAGLib diagonalization */ 
{ 

f ree_dvector (mat_l , ,diiiil*diml-l) ; 
f re6_dvector Ce_l ,0 ,diml-l) ; 
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free_dvector(r_l,0,diml-l) ; 

f ree_dvector (v_l , , diml*diml-l) ; 

f ree_dvector (mat_2 , , dim2*dim2-l) ; 

f ree_dvector (e_2,0,dim2-l) ; 

f ree_dvector (r_2,0,diin2-l) ; 

f ree_dvector (v_2 , ,dim2*dim2-l) ; 

if (dim3>0) 
{ 

f ree_dvector (mat_3 , , dlm3*dim3-l) ; 

f ree_dvector (e_3,0,dim3-l) ; 

f ree_dvector (r_3,0,diin3-l) ; 

f ree_dvector ( v_3 , , diin3*dim3- 1 ) ; 

f ree_dvector (mat_4 , , dim4*dim4-l ) ; 
free_dvector(e_4,0,dim4-l) ; 
free.dvector (r_4,0,dim4-l) ; 
f ree_dvector (v_4 , , dim4*dim4-l) ; 

} 
} 

void free_ma(int dim) 

/* deallocate matrices used for asymmetric NAGLib diagonalization */ 
{ 

free_dvector(nia,0,dim*dim-l) ; 
free_dvector(e,0,di!ii-l) ; 
f ree_dvector (r ,0,diin-l) ; 
f ree_dvector (v,0,dim*dim-l) ; 

} 

void read_input() 

/* read in inputfile created by 'meta_mesonix. c ' */ 
{ 

FILE *fp; 

char *name_IN = "Input_mes.dat"; 

char comment [80]; /* dummy variable to read out comments*/ 



fp = fopen(name_IN,"r"); /* read in INPUT_MES.DAT */ 



fscanf (fp,"'/.s" 


, comment) ; 












fscanf (fp,"'/.s' 


, comment) ; 












fscanf (fp,"y.s' 


, comment) ; 












fscanf (fp,"y.s' 


, comment) ; 












fscanf (fp,"7.1f", fern) ; 


/* 


read 


in 


<fermion mass> */ 




fscanf (f p, "7,s' 


, comment) ; 












fscanf (fp, ""/.If 


",&lambda) ; 


/* 


read 


in 


<cut-off> */ 




fscanf (fp/'/'.s' 


, comment) ; 












fscanf (fp, '"/.If 


",&ALPHA) ; 


/* 


read 


in 


<coupling> */ 




fscanf (f p, "y.s' 


, comment) ; 












fscanf (f p, "y.d' 


,&test) ; 


/* 


read 


in 


<test variable> */ 




fscanf (f p, "y.s' 


, comment) ; 












fscanf (f p, "y.d' 


,&J_z) ; 


/* 


read 


in 


<J_z> */ 




fscanf (f p, "'/.s' 


, comment) ; 












fscanf (f p, "y.d' 


,&Anni) ; 


/* 


read 


in 


control for annih. 


graph */ 


fscanf (fp,"y.s' 


, comment) ; 












fscanf (fp,"y.d' 


,&N1); 


/* 


read 


in 


# mesh points <N1> 


*/ 


fscanf (fp, "y,s' 


, comment) ; 












fscanf (fp, "y.d' 


,&N2); 


/* 


read 


in 


# mesh points <N2> 


*/ 


fscanf (fp, "y.s' 


, comment) ; 












fscanf (fp,"y.d' 


, <mumber_EV) ; 


/* 


read 


in 


# of EVs written in file */ 


fscanf (fp,"y.s" 


, comment) ; 












fscanf (fp, "y.d' 


,&trix) ; 


/* 


read 


in 


control for Ctrick 


ab initio*/ 


fscanf (f p, "y.s' 


, comment) ; 












fscanf (f p, "'/.s' 


,def ault_directory) ; 


/* read 


in <default directory> */ 


fscanf (fp, "y.s' 


, comment) ; 
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f scanf (fp, "%s" ,special_id) ; /* read in <special identif ication> */ 
f scanf (fp, "/is" .comment) ; 

f scanf (fp, "/is" ,&print_EF) ; /* read in control for writing */ 

/* eigenfctns in file */ 

f close (f p) ; 

name.EV = (char *)malloc(120) ; 

sprintf (name_EV, "'/.sEigenValues/EVs'/.s J'/.d . dat " , 

default_directory,special_id, J_z) ; /* create filename for EVs */ 



void logixO 
/* write LOG 



file */ 



name_log = (char *)malloc(120) ; /* name for LOG file 
sprintf (naine_log, "'/.smesonix . log" .default .directory) ; 
f-log = f open(name_log, "a") ; 



*/ 



fprintf (f _log, 
fprintf (f _log, 
fprintf (f.log, 
fprintf (f.log, 
fprintf (f.log, 
fprintf (f.log, 
fprintf (f _log, 
fprintf (f _log, 
fprintf (f _log, 
fprintf (f _log, 
fprintf (f _log, 
fprintf (f _log, 
fprintf (f _log, 
fprintf (f _log, 
f close (f.log) ; 



\nM E S N I X 



PARAMETERS :\n") ; 



m 




%18 


lambda 




%18 


alpha 




•/.18 


:\n"); 






test 




y.5d 


J_z 




'/.5d 


Anni 




'/.5d 


Nl 




'/.5d 


N2 




'/.5d 


number.EV 




•/.5d 



LOG file\n"); 
\n\n"); 



12f ; f ermion mass\n" ,m) ; 
12f ; cut-of f\n", lambda) ; 
12f ; coupling\n", ALPHA); 

test variable\n" ,test) ; 
ang. momentum sector \n" , J.z) ; 
annihilation graph\n" , Anni) ; 
number of MU values\n" ,N1) ; 
number of THETA values\n" ,N2) ; 
number of eigenvalues\n" , number.EV) ; 



\n- 



-\n\n"); 



void close_log_f ile 

/* final comment and close of LOG file */ 



fprintf (f.log, "\n \n\n") ; 

f close (f.log) ; 



void allocationO 

/* allocate matrices used */ 

{ 

init_ma (4*N1*N2) ; /* initialize matrix for NAGLIB diag */ 

if ((testai6384) !=0) init_mat (N_sym_l ,N_sym_2 , ,0) ; 
else init.mat(Nl*N2,Nl*N2+Nl,Nl*N2,Nl*N2-Nl) ; 

/* initialize matrices 1-4 for NAGLIB diag. */ 
mu = dvector (1 ,NN1+1) ; 

theta = dvector (1,NN2+1) ; 
wmu = dvector(l,NNl+l) ; 
wtheta = dvector(l,HN2+l) ; 

} 



void deallocationO 

/* deallocate matrices used */ 

{ 

if ((test&16384) !=0) free.mat(N.sym.l,M.sym.2,0,0) ; 
else free.mat(Nl*N2,Nl*N2+Nl,Nl*H2,Nl*N2-Nl) ; 
f ree.ma (4*N1*N2) ; 

} 



void mainO 
{ 

int i , j ; 
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FILE *fp; 

pr intf ( " \nMESONIX \n" ) ; 

read_input() ; /* read in INPUT_MES.DAT */ 

loglxO ; /* write header of LOG file */ 

allocationO; /* allocate all objects used */ 

p_bohr = m* ALPHA/2; /* Bohr momentum */ 

if (old_CT==0 ) coulomb_trix(trix,3,Nl) ; /* init Coulomb trick */ 

if (trix==0) /* if not C. trick from scratch, then . . .*/ 

{ 

fp = f open(name_EV, "a") ; /* open file for eigenvalues */ 
fprintf (fp,"'/.3d\n",Nl) ; 

f close(fp) ; 

f_log = f open(naine_log, "a") ; /* open LOG file */ 
fprintf (f .log, "Time for\n"); 
f close (f_log) ; 

numerix(Nl,N2) ; /* calculate weights for numerical integration */ 

f-log = f open(name_log, "a") ; /* open LOG file for subroutines */ 
plot_running_coupling() ; 

if ((tests (4096+16384)) == 0) /* J_z=0 */ 
{ 

genetixO; /* set up Hamilton matrix */ 

arithmetix(Nl*N2,Nl*N2+Nl,Nl*N2,Nl*N2-Hl) ; /* diagonalize */ 

} 

else /* J.zOO */ 
{ 

if ((test&16384) == 0) /* standard, unsynmietrized Hamiltonian */ 
{ 

asymmetrixO ; /* set up Hamiltonian without symmetries */ 
arithmetix(4*Nl*N2, 0,0,0) ; 

} 

} 

publixO; /* write eigenvalues in file */ 

} 

close_log_f ileO ; 

deallocationO ; /* deallocate all objects used */ 

printf ("\n \n"); 
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Numerix 



/* */ 

/* program :NUMERIX */ 

/* */ 

/* class : subroutine of MESONIX */ 

/* */ 

/* file name : "numerix. c" */ 

/* */ 

/* purpose : NUMERIX calculates the weights and abscissae */ 

/* for numerical integration according to the */ 

/* Gauss-Legendre-algorithm. */ 

/* */ 

/* */ 

/* */ 

/* Subroutines : dOlbcf (NAGLib) for Gauss-Legendre discretization */ 

/* */ 

/* Global variables : mu[l..NNl] = abscissa values in mu direction */ 

/* wmu[l..NNl] = weights in mu direction */ 

/* thetaCl . .NN2] = abscissa values in cos(theta) */ 

/* direction */ 

/* wmu[l..NN2] = weights in cos(theta) direction */ 

/* */ 



extern double d01bcf_( int *itype, double *a, double *b, double *c, double *d, 
int *n, double *wmu, double *mu, int *if ail) ; 

/* NAGLib: Gauss-Legendre discretization */ 

void numerixC int nl, int n2 ) /* nl = number of mu-values */ 

/* n2 = number of theta-values */ 
/* test == 1 prints values in FILE */ 

/* main function: calculate weights and abscissae */ 

{ 

int i ; 

int itype.ifail; /* for dOlbcf */ 

double a,b,c,d; /* for dOlbcf */ 

double *weight ,*abscis; /* for dOlbcf */ 
char *name; /* name for file in case test=l */ 
FILE *fp; 

timeC&timel) ; 



printfC > NUMERIXC Ml='/.2d, N2 ="/.2d)\n", 

nl,n2) ; 

f-log = f open(name_log, "a") ; 

/* +++++ Discretize in mu-direction, boundaries [0, lambda/2] ++++ */ 

weight = dvector (0,nl-l) ; /* initialize fields for dOlbcf */ 

abscis = dvector (0,nl-l) ; 

itype =0; /* parameters for dOlbcf */ 



a = l/(l+lambda/2.0/p_bohr) ; 
b = c = d = 1.0; 

dOlbcf _(&itype, &a, &b, &c, &d, &nl, weight, abscis, &if ail) ; 

/* Gauss-Legendre discretization */ 

for (i=l; i<=nl; i++) 
{ 

wmu[i] = weight [i-1] /abscis [i-1] /abscis [i-l]*p_bohr; 
mu[i] = (l/abscis[i-l]-l)*p_bohr; 

} 

free_dvector(weight,0,nl-l) ; 
free_dvector(abscis,0,nl-l) ; 
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/* ===> Discretize in cos(theta)-direction, boundaries [-1,1] <=== */ 



weight = dvector (0,n2-l) ; /* initialize fields for dOlbcf */ 

abscis = dvector (0,n2-l) ; 

itype =0; /* parameters for dOlbcf */ 

a = -1.0; 

b = c = d = 1.0; 



dOlbcf _(&itype, &a, &b, &c, &d, &n2, weight, abscis, &if ail) ; 

/* Gauss-Legendre discretization */ 

for (i=l; i<=n2; i++) 
{ 

theta[i] = abscis [i-1]; 
wtheta[i] = weight [i-1]; 

} 

free.dvector (weight, 0,n2-l) ; 
free.dvector (abscis, 0,n2-l) ; 

/* =========> test==l : write Gauss-Legendre values in a file <======= */ 

if ((test&l) != ) 
{ 

name = (char *)inalloc(120) ; 

sprintf (name , "'/,sgauss_values . dat " ,def ault_directory) ; /* filename */ 
fp = f open(naine, "w") ; 

fprintf (fp, "Gauss-Legendre integration : MESONIX.C\n") ; 

fprintf (fp," \n\n") ; 

fprintf (fp," mu- integration 

theta-integration\n\n") ; 
fprintf (fp," n abscissa weight 

abscissa weight\n"); 

fprintf (fp, " 

\n"); 

for (i=l; i<=nl; i++) 
{ 

fprintf (fp,"/'.3d I •/.12.8f •/.12.8f I •/.12.8f 

■/.12.8f \n",i,mu[i] ,wmu[i] ,theta[i] ,wtheta[i] ) ; 

} 

fprintf (fp," 

\n"); 

f close (fp) ; 

} 

/* =======> write time used for calculation in LOG file <========== */ 



time (&time2) ; 

fprintf (f_log, " ...weights: '/,d sec\n" ,time2-timel) ; 

f close (f _log) ; 
return ; 
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Coulomb _trix 

/* */ 

/* Program : COULOMB.TRIX */ 

/* */ 

/* Class : Subroutine */ 

/* */ 

/* Purpose : COULOMB.TRIX creates the files for MESONIX */ 

/* which contain the values for the different */ 

/* Coulomb tricks, for general J_z. */ 

/* */ 

/* */ 

/* */ 
/* Names of the generated files: */ 
/* */ 
/* 1) CTrick_J<n>_uu.dat \ */ 
/* 2) CTrick_J<n>_ud.dat \_ n = ...-2,-1,0,1,2,... */ 
/* 3) CTrick_J<n>_du.dat / */ 
/* 4) CTrick_J<n>_dd.dat / */ 
/* */ 


/* External functions: */ 


/* */ 

/* dOlalf NAGLib integration of function with singularities */ 

/* dOlajf NAGLib integration of function without sing. */ 

/* */ 

extern double d01alf_( double (*funktion) (double *) , double *a, double *b. 



extern double d01ajf_( double (*funktion) (double *) , double *a, double *b, 
double *epsabs, double *epsrel, double *result, 
double *abserr, double w[lp], int *lm, int iw[lip]. 



int *npts, double points [np], double *epsabs , 
double *epsrel, double *result, double *abserr, 
double w[lp],int *lm, int iw[lip] , int *liw, 
int *if ail) ; 
/* integration routine from NAGLib */ 

extern double d01ajf_( double (*funktion) (double *) , double *a, double 
double *epsabs, double *epsrel, double *result, 
double *abserr, double w[lp], int *lm, int iw[l 
int *liw, int *if ail) ; 

/* integration routine from NAGLib */ 

double old_continous(int jm, int jt) 
/* old Krautgaertner Coulomb trick */ 
{ 

int ifail,npts,iw[lip] ,lw,liw; 

double a, b, points [2] ,w[lp] ,abserr,epsabs,epsrel,coul_conti; 

a = 0.0; 

b = lambda/2.0; 

jjl = jm; 
jj2 = jt; 

/* NAGLIB-ROUTINE */ 

epsabs = le-6; 
epsrel = le-8; 
npts = 1; 

points [0] = mu[jm]; 
liw = lip; 
Iw = Ip; 
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if ail = -1; 

d01alf_ (integrand, Sia,&b ,&npts .points ,&epsabs ,&epsrel, 

&coul_cont i , feabserr , w , &lw , iw , &liw , &if ail) ; 
return coul_conti; 



double Coulomb_trlck_f unction (double *theta_lokal) 

/* continous C. trick coimterterm, to be integrated over <theta> */ 

/* for fixed <mu> */ 

{ 

/* Parameters = global variables are: 
mu_global 

Jjl.jj2 
J_z_Coulomb 
spins_parallel */ 
double result, thetaa; 

thetaa = *theta_lokal; 
result =0.0; 

x[l] = (1.0 + mu_global*thetaa/sqrt (m*m + mu_global*mu_global))/2.0; 
k[l] = mu_global*sqrt (1 .0-thetaa*thetaa) ; 

x[2] = (1.0 + mu[jjl]*theta[jj2]/sqrt(m*m + mu[j j 1] *mu[j j 1] ) )/2 . 0; 
k[2] = mu[jjl]*sqrt(1.0-theta[jj2]*theta[jj2]); 

if (( fabs(x[l]-x[2]) >= le-8 I I f abs (k [1] -k [2] ) >= le-8 ) 
&fe (fabs(k[l]) >= le-8) && (fabs(k[2]) >= le-8)) 

{ 

aa = (x[l]-x[2])*(x[l]-x[2])*m*ni/2.0*(l/(l-x[2])/(l-x[l]) 
+l/x[i]/x[2]) + k[l]*k[l] + k[2]*k[2] 
+ (x[l]-x[2])/2*(k[l]*k[l]*(l/(l-x[l])-l/x[l]) - 
k [2] *k [2] * (1/ (1-x [2] ) -1/x [2] ) ) ; 

if ( (aa*aa-4*k [1] *k [1] *k [2] *k [2] ) >0 . 0) 
{ 

A = 1.0/sqrt( aa*aa - 4*k[l]*k[l]*k[2]*k[2] ); 
B = (1.0-aa*A)/2.0; 

if (spins_parallel==l) result = Gl (2 , 1 , J_z_Coulomb) ; 
else result = G2(2,l,J_z_Coulomb) ; 

result *= 2.0*mu_global*mu_global*x [1] *(1 .0-x [1] )/ 
sqrt(in*in + inu_global*mu_global) * 

( 1+mu [ j j 1] *inu [ j j 1] *inu [ j j 1] /p_bohr/p_bohr/ p_bohr/8 . 0) 

/ (l+mu_global*mu_global*mu_global 
/p_bohr/p_bohr/p_bolir/8 . 0) /PI ; 
} 

} 

return result; 
} 

double Coulomb_trick_integrand (double *mu_lokal) 

/* C. trick counterterm, already integrated over <theta>, dependent on <mu> */ 
{ 

/* Parameters = global variables are 
mu.global 

J_z_Coulomb 
spins_parallel */ 

int ifail,npts,iw[lip] ,lw,llw; 

double a, b, points [2] ,w[lp] ,abserr,epsabs,epsrel,result_coul; 
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mu_global = *mu_lokal; 
a =-1.0; 
b = 1.0; 
epsabs = le-12; 
epsrel = ie-3; 
npts = 1; 

points [0] = theta[jj2] ; 
liw =lip; 
Iw = Ip; 
if ail =-1; 

if (f abs (mu_global - mu[jjl]) <= le-10) 
{ 

dOlalf _ (Coulomb_trick_f unction , &a , &b , &npts , points , feepsabs , 
feepsrel , &result_coul , ftabserr , w , &lw , iw , ftliw , &if ail) ; 

} 

else 

dOlaj f _ (Coulomb_trick_f unction , &a , &b , feepsabs , feepsrel , 
&result_coul , feabserr , w , felw, iw, feliw , feif ail) ; 
return result.coul; 



double Coulomb_discrete(int jm, int si, int s2) 

/* discrete part of Coulomb trick, called by PHYSIX */ 

{ 

int km , kt ; 

double result , spinf unct ion ; 



result =0.0; 

for ( km = 1; km <= Nl; ++km ) 
{ 

for ( kt = 1; kt <= N2; ++kt ) 
{ 

x[l] = (1.0 + mu[km]*theta[kt]/sqrt(m*m + mu[km] *mu[km] ) )/2 . 0; 
k[l] = mu[km]*sqrt(1.0-theta[kt]*theta[kt]); 

if ( fabs(x[l]-x[2]) >= le-8 II f abs(k[l] -k[2] ) >= le-8 ) 
{ 

aa = (x[l]-x[2])*(x[l]-x[2])*m*m/2.0*(l/(l-x[2])/(l-x[l]) 
+l/x[l]/x[2]) + k[l]*k[l] + k[2]*k[2] 
+ (x[l]-x[2])/2*(k[l]*k[l]*(l/(l-x[l])-l/x[l]) - 
k [2] *k [2] * (1/ (1-x [2] ) -1/x [2] ) ) ; 

A = 1.0/sqrt( aa*aa - 4*k[l]*k[l]*k[2]*k[2] ); 

B = (1.0-aa*A)/2.0; 



if (sl>0) 
{ 

if (sl==s2) spinfunction = -Gl(2,l,J_z); 

else spinf imction = -G2(2, 1 , J_z) ; 

} 

else 
{ 

if (sl==s2) spinfunction = -Gl (2, 1 ,-J_z) ; 

else spinfunction = -G2(2, 1 ,-J_z) ; 

} 

if (old.CT !=0) spinfunction = -G2(2, 1 , J_z) ; 

result += 1.0/PI*wtheta[kt] *wmu [km] *spinf unction 
*2 . 0*mu [km] *mu [km] *x [1] * (1 . 0-x [1] ) 
/sqrt(m*m + mu [km] *mu [km] ) * 

( 1+mu [ jm] *mu [ jm] *mu [ jm] /p_bohr/p_bohr/p_bohr/ 
/ ( 1+mu [km] *mu [km] *mu [km] /p_bohr/p_bohr /p_bohr/ 8) ; 



} 

} 

} 

return result; 
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} 

void install_Coulomb_trick() 

/* read data from 'CTrick_Jn_ss.dat' into CT[] */ 
{ 

int i,j ,k,l,lo[2] ,hi[2] .error; 
double value; 
FILE *f p [4] ; 

error = 0; 

for (i=0; i<=3; ++i) 
{ 

fp[i] = f open(CT_naine [i] , "r") ; 
fscanf (fp[i] , "y.d" ,&lo [0] ) ; 
fscanf (fp[i] ,"y.d",&lii [0] ) ; 

if (i>0 && (lo[l]!=lo[0] II hi[l] !=hi[0])) error = 1; 

lo[l]=lo[0]; 

hi[l]=hi [0]; 

} 

if (error ! =0) 

printf("»»> ERROR(install_Coulomb_trick) : lo[i] !=lo[j]\n") ; 
if (hi[l]<Nl) 

printf("»»> ERROR(install_Coulomb_trick) : hi[l] <Nl\n"); 

for (i=lo[l]; i<Nl; i+=2) /* ignore first data until actual Nl */ 
{ 

for (1=0; 1<=3; ++1) f scanf (f p [1] , '"/.If " ,&k) ; 

for (j=l; j<=i*i; 

{ 

for (1=0; 1<=3; ++1) fscanf (fp[l] , '"/.If", fevalue) ; 

} 

} 

for (1=0; 1<=3; ++1) f scanf (fp[l] , '"/.If " ,&k) ; 

for (i=0; 1<=N1-1; /* read out needed data */ 

for (j=0; j<=N2-l; 

{ 

for (1=0; 1<=3; ++1) 
{ 

fscanf (fp[l] .'"/.If ".fevalue) ; 
CT[1] [i] [j] = value; 

} 

} 

for (1=0; 1<=3; ++1) f close (fp[l] ) ; 

} 

void special_install_Coulomb_trick() 

/* read data from 'CTrick_Jn_ss.dat' into CT[] , special case */ 
{ 

int i, j ,k,l,nl [2] ,n2[2] .error; 
double value ; 
char * comment [80] ; 
FILE *fp[4]; 

error = 0; 

for (i=0; i<=3; ++i) 
{ 

fp[i] = fopen(CT_name[i],"r"); 
fscanf (fp[i] . ""/.s" .comment) ; 
fscanf (fp [i] . "%d" .ftnl [0] ) ; 
fscanf (fp [i] . '"/.d" .&n2 [0] ) ; 

if (i>0 «E& (nl[l] !=nl[0] II n2[l] !=n2[0] )) error = 1; 

nl[l]=nl [0] ; 
n2[l]=n2[0] ; 

} 

if (error !=0) 
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printf("»»> ERROR(install_Coulomb_trick) : lo [i] !=lo [j] \n") ; 
if (nl[l]<Nl) 

printf("»»> ERROR(install_Coulomb_trick) : nl[l] < Nl\n"); 

for (i=0; i<=Nl-l; ++i) /* read out needed data */ 
for (j=0; j<=N2-l; ++j) 
{ 

for (1=0; 1<=3; ++1) 
{ 

fscanf (fp[l] ,"y.lf ".fevalue) ; 
CT[1] [i] [j] = value; 

} 

} 

for (1=0; 1<=3; ++1) f close (fp[l] ) ; 

} 



void special_coulomb_trix(int f rom.scratch) 
/* calculate or initialize the Coulomb coimterterms */ 
/* !!!!!!!!!!! special case Nl ! =N2 !!!!!!!!!!!!! */ 


/* from_scratch != — > really calculate the terms */ 
/* == — > just install the trick */ 

{ 

int i, j ,1,NNN, counter; 

int if ail ,npts , iw[lip] , Iw, liw; 

double a, b, points [2] ,w[lp] ,abserr ,epsabs,epsrel , coul_result ; 
FILE *f p [4] ; 

printf("\n (special)\n") ; 

for (i=0; i<=3; ++i) /* create file names */ 
{ 

CT_name[i]=(char *)malloc(80); 

sprintf (CT.name [i] , "•/.sCoulomb.Trick/sCTrickZsJ" , 

default_directory,special_id) ; 
switch (i) 
{ 

case : sprintf (CT.name [i] , " '/,s'/.d_uu . dat " , CT_name [i] , J_z) ; 
break; 

case 1 : sprintf (CT.name [i] , ""/.s'/.d.ud.dat" ,CT_name[i] , J_z) ; 
break; 

case 2 : sprintf (CT.name [i] , " '/.s'/.d.du . dat " , CT.name [i] , J_z) ; 
break; 

case 3 : sprintf (CT.name [i] , ""/.s'/.d.dd.dat" , CT.name [i] , J.z) ; 
break; 

} 

} 

if (f rom.scratch == 0) special.install.Coulomb.trickO ; 

else 
{ 

for (i=0; i<=3; ++i) /* open file */ 
{ 

printf("»»> INFO: Opening file ''/.s' \n" .CT.name [i] ) ; 
fp[i] = fopen(CT.name[i] , "w") ; 

} 

for (i=0; i<=3; ++i) fprintf (fp [i] , "Special : \ny,2d\ny,2d\n" ,N1,N2) ; 

numerix(Nl,N2) ; /* calculate Gauss.Legendre values */ 

for (i=l; i<=Nl; ++i) 
{ 

printfC i = '/.2d\n",i); 

for (j=l; j<=N2; 

{ 
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printfC j = •/.2d\n",j); 

covinter = 0; 

for (1=1; 1>=-1; l-=2) 

for (spins _parallel=l ; spins _parallel>=-l ; 

spins _parallel-=2) 

{ 

J_z_Coulomb = 1*3 _z; 

x[2] = (1.0+mu[i]*theta[j]/sqrt(in*m+mu[i]*mu[i]))/2.0; 

k[2] = mu[i]*sqrt(1.0-theta[j]*theta[j]) ; 

jjl = i; 

jj2 = j; 

a = 0.0; 

b = lambda/2.0; 

epsabs = le-6; 

epsrel = le-8; 

npts = 1; 

liw = lip; 

Iw = Ip; 

if ail = -1; 

points [0] = mu[i]; 

if ( J_z_Coulomb ! =0 I I 1>0) 

{ 

dOlalf _(Coulomb_trick_integrand,&a,&b,topts, 
points , &epsabs , ftepsrel , &coul_result , 
&abserr,M,&lM,iw,&liM,&ifail) ; 
f printf (f p [counter+spins_parallel* (1-1) /2] , 

"•/.18 . 12f \n" , coul_result) ; 
if (J_z_Coulomb==0) 

f printf (f p [coimter+2+spins_parallel] , 
""/.18.12f\n",coul_result); 

} 

++coimter ; 

} 

} /* j loop */ 
} /* i loop */ 

for (i=0; i<=3; ++i) f close (fp [i] ) ; 
} /* END: from.scratch != */ 

} 



void coulomb_trix(int f rom_scratch, int lo_N, int hi_N) 
/* calculate or initialize the Coulomb counterterms */ 
/* */ 

/* f roni_scratch != — > really calculate the terms */ 
/* == — > just install the trick */ 

/* lo_N — > lowest discretization */ 
/* hi_N — > highest discretization */ 
{ 

int i, j ,l,n,NNN, counter; 

int if ail, npts, iw [lip] ,lw, liw; 

double a, b, points [2] ,w[lp] ,abserr, epsabs, epsrel, coul_result; 
FILE *f p [4] ; 

printf ("\n > COULOMB_TRIX\n") ; 

if (N1!=N2) special_coulomb_trix(f rom_scratch) ; 

else 

{ 

for (i=0; i<=3; ++i) /* create file names */ 
{ 

CT_name[i]=(char *)malloc(80) ; 

sprintf (CT.name [i] , "'/.sCoulomb.Trick/CTrick'/.s J" , 

default .directory, special.id) ; 
switch (i) 
{ 
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case : sprintf (CT.name [i] , " '/.s'/.d.uu . dat " , CT_name [i] , J_z) ; 
break; 

case 1 : sprintf (CT.name [i] , ""/.s'/.d.ud.dat" ,CT_name[i] , J_z) ; 
break; 

case 2 : sprintf (CT.name [i] , " '/.s'/.d_du . dat " , CT.name [i] , J_z) ; 

break; 

case 3 : sprintf (CT.name [i] , ""/.s'/.d.dd.dat" , CT.name [i] , J.z) ; 
break; 

} 

> 

if (f rom.scratch == 0) install. Coulomb.trickO ; 

else 
{ 

for (i=0; i<=3; ++i) /* open file */ 
{ 

printf("»»> INFO: Opening file ''/.s' \n" .CT.name [i] ) ; 
fp[i] = fopen(CT.name[i] ,"w"); 

} 

for (i=0; i<=3; ++i) f printf (f p[i] , "•/,2d\n'/.2d\n" ,lo.N,hi.N) ; 

for (n=lo.N; n<=hi.N; n+=2) 

{ 

numerixCn.n) ; /* calculate Gauss. Legendre values */ 

printf ("Counterterms: n/hl_N = •/.2d/y,2d\n" ,n,hi_N) ; 
for (i=0; i<=3; ++i) fprintf (f p [i] , "•/.2d\n" ,n) ; 
for (i=l; i<=n; ++i) 
{ 

printf (" i = •/.2d\n",i); 

for (j=l; j<=n; 

{ 

printf (" j = •/.2d\n",j); 

counter = 0; 

for (1=1; 1>=-1; l-=2) 

for (spins .parallel=l; spins.parallel>=-l;spins.parallel-=2) 
{ 

J.z.Coulomb = l*J.z; 

x[2] = (1.0+mu[i]*theta[j]/sqrt(m*m+mu[i]*mu[i]))/2.0; 
k[2] = mu[i]*sqrt(1.0-theta[j]*theta[j]) ; 

jjl = i; 

jj2 = j; 

a = 0.0; 

b = lambda/2.0; 

epsabs = le-6; 

epsrel = le-8; 

npts = 1; 

liw = lip; 

Iw = Ip > 

if ail = -1; 

points [0] = mu[i]; 

if (J_z_Coulomb !=0 I I 1>0) 

{ 

dOlalf . (Coulomb.trick_ integrand , &a , &b , &npt s , 
points , ftepsabs , feepsrel , &coul_result , 
&abserr,w,&lw,iw,&liw,&ifail) ; 

fprintf (f p [counter+spins .parallel* (1-1) /2] , 
"•/.18.12f\n",coul.result) ; 

if (J.z.Coulomb==0) 

fprintf (f p [counter+2+spins.parallel] , 
"•/.18 . 12f \n" , coul.result) ; 

} 

++counter; 

} 

} 

} 

} 
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for (i=0; i<=3; f close (fp [i] ) ; 

} /* END: from.scratch != */ 
} /* END: N1==N2 */ 

} 



double Integrand ( double *mu_c ) 

/* Krautgaertner counterterm, analytically integrated over <theta> */ 

/* (only used if <old_CT>!=0) */ 

{ 

double el,elp,ell,f,g,h,kl,a,b,c , d,x_cc,k_cc,mu_cc ; 
double xpl,xp2,xp3,xml,xm2,xm3,abk; 
double ql,q2,q3,q4,q5; 

mu_cc = *mu_c; 

x_cc = (1.0 + mu[jjl]*theta[jj2]/sqrt(m*m + mu[j j 1] *mu[j j 1] ))/2.0; 
k_cc = mu[j jl]*sqrt(1.0 - theta[j j2] *theta[j j2] ) ; 

el = sqrt(m*m + mu[j jl] *mu[j jl] ) ; 
elp = sqrt(m*m + mu_cc*mu_cc) ; 

kl = elp/el +el/elp; 

a = mu[j jl] *mu[j jl] *mu_cc*mu_cc*(4. +theta[j j2] *theta[j j2] * 
(elp - el)*(elp-el)/el/elp/el/elp) ; 

b = -2 . 0* (mu [ j j 1] *mu [ j j 1] + mu_cc*mu_cc)*mu[j j 1] *mu_cc*theta[j j2] *kl ; 

c = (inu[j jl] *mu[j jl] + inu_cc*mu_cc) * (mu [j j 1] *mu[j j 1] + inu_cc*mu_cc) 
-4.0*mu[j jl] *mu[j jl] *mu_cc*inu_cc*(l .0-theta[j j2] *theta[j j2] ) ; 

d = m*iii/(1.0 - x_cc)/x_cc/elp + 2.0*mu_cc*mu_cc/elp 
+ k_cc*k_cc/2.0/elp/x_cc/(1.0 - x_cc) 

+ 0.5/(1.0 - x_cc) /x_cc/elp*(mu[j j 1] *mu[j j 1] + mu_cc*mu_cc) ; 

ell = 2.0*m*m/(1.0 - x_cc)/x_cc/elp/elp*mu_cc*(x_cc - 0.5) 
1.0/x_cc/(1.0 - x_cc)/elp* 
(-mu [ j j 1] *mu_cc*theta [ j j 2] *kl*0 . 5 + 

(mu [ j j 1] *mu [ j j 1] + mu_cc*mu_cc)*mu_cc/elp*(x_cc - 0.5)); 

f = -2.0*mu_cc*mu_cc/elp - 0.5*mu_cc*mu_cc/x_cc/(1.0 - x_cc)/ 

elp/elp/elp*k_cc*k_cc - 1 . 0/elp/elp/x_cc/(l . - x_cc)* 
mu_cc*(x_cc - 0. 5) *mu [j j 1] *mu_cc*theta[j j2] *kl ; 

h = -0.5/x_cc/(1.0 - x_cc)/elp; 

xpl = a + b + c; 

xp2 = (mu_cc*mu [j j 1] *theta[j j2] *kl - (mu [j j 1] *mu [j j 1] + mu_cc*mu_cc) ) * 
(mu_cc*mu[j j 1] *theta[j j2] *kl - (mu[j j 1] *mu[j j 1] + mu_cc*mu_cc) ) ; 
xp3 = fabs(mu_cc*mu[jjl]*theta[jj2]*kl-(mu[jjl]*mu[jjl] + mu_cc*mu_cc) ) ; 

xml = a - b + c; 

xm2 = (mu_cc*mu[j jl] *theta[j j2] *kl + (mu [j j 1] *mu [j j 1] + mu_cc*mu_cc) ) * 
(mu_cc*mu[j jl] *theta[j j2] *kl + (mu [j j 1] *mu [j j 1] + mu_cc*mu_cc) ) ; 

xm3 = fabs(mu_cc*mu[jjl]*theta[jj2]*kl + (mu [ j j 1] *mu [ j j 1] 
+ mu_cc*mu_cc) ) ; 

if ( (2.0*sqrt(a)*xm3 - 2.0*a + b == 0) 

I I (2.0*sqrt(a)*xp3 + 2.0*a + b == 0)) return 0; 

else 
{ 

abk = 1.0/sqrt(a)*log(fabs((2.0*sqrt(a)*xp3 + 2.0*a + b) 

/(2.0*sqrt(a)*xm3 - 2.0*a + b))); 
return -1.0/PI* 
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(d*abk*mu_cc*mu_cc + ell*mu_cc*mu_cc/a 
*(sqrt(xp2) - sqrt(xm2) - 0.5*b*abk) 

+f * ( (mu_cc*mu_cc*0 . 5/a - 0.75*b*mu_cc*mu_cc/a/a) *sqrt (xp2) 
-(-mu_cc*mu_cc*0.5/a - 0.75*b*mu_cc*mu_cc/a/a)*sqrt(xm2) 
+(3.0*b*b - 4.0*a*c)*abk/8.0*mu_cc*mu_cc/a/a) 

+ 2.0 *h*mu_cc*mu_cc) 

* (1 . 0+mu [j j 1] *mu [j j 1] *mu[j j 1] /p_bohr/p_bohr/p_bohr/8 . 0) 
/(I .0+mu_cc*mu_cc*mu_cc/p_bohr/p_bohr/p_bohr/8. 0) ; 

} 

} 

Genetix 

/* */ 

/* program :GENETIX */ 
/* */ 
/* class : subroutine of MESONIX */ 
/* */ 
/* file name : "genetix. c" */ 
/* */ 
/* purpose : GENETIX generates the Hamilton matrix of the */ 
/* system. It uses the subroutine PHYSIX for */ 
/* calculation of the single elements. */ 
/* */ 


/* */ 

/* External functions: */ 

/* */ 

/* */ 

/* physix subroutine of 'mesonix', calculates one matrix */ 

/ * element */ 

/* */ 

void genetixO 

/* create 4 Hamiltonians, symmetries C,H */ 
{ 



double f actor_jl,f actor_j2,factor_j3; /* factors for wavefunctions */ 
double f actor_il , f actor_i2 , f actor_i3 ; 

int jnum_l, jnum_2, jnum_3, jnum_4; /* actual Hamiltonian row */ 

int inum_l , inum_2, inum_3, inum_4; /* actual Hamiltonian column */ 

time(&timel) ; 

printfC > GENETIX \n"); 

Sqrt2 = l/sqrt(2.0); 

jnum_l = jnum_2 = jnum_3 = jnum_4 = 0; 
sjl = sil = 1; 

/* =====================> rows of Hamiltonian <================= */ 

for ( jm=l; jm <= Nl; ++jm) /* loop for <mu> */ 
{ 

for ( jt=l; jt <= (N2+l)/2; ++jt) /* loop for <theta> */ 
{ 
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printf ("jm/Nl =%2d/"/.2d jt/N2 ="/.2d/y,2(i\n" , jm.Hl , jt , (N2+l)/2) ; 
for ( sj2 = -1; sj2 <= 1; sj2 += 2 ) /* loop for <spiii> */ 
{ 

factor_jl = factor_j2 = factor_j3 = 1; 

if ( 2*jt-l != N2 ) 
{ 

++jnum_l ; 
++jnuin_2; 
++jiiujn_3; 
++jnum_4; 



j .parity = 0; 

} 

else 
{ 

if (sjl == sj2) 
{ 



/* not end of loop */ 



} 

else 
{ 



++jnimi_l; 
++jnimi_2; 

factor_jl = factor_j2 = Sqrt2; 

j.parity = 1; /* end of loop && sl=s2 */ 



++jnum_2 ; 
++jnum_3 ; 

factor_j2 = factor_j3 = Sqrt2; 

j_parity = -1; /* end of loop && sl!=s2 */ 



/* ===============> columns of Hamiltonian <= 

inum_l = inum_2 = inum_3 = inum_4 = 0: 



*/ 



for ( im=l; im <= Nl; ++im) /* loop for <mu> */ 
{ 

for ( it=l; it <= (N2+l)/2; ++it) /* loop for <theta> */ 
{ 

for ( si2 = -1; si2 <= 1; si2 += 2 ) /* loop for <spin> */ 
{ 

factor_il = factor_i2 = factor_i3 = 1; 



if ( (2*it-l) 
{ 

++inum_l ; 
++inum_2 ; 
++inum_3 ; 
++inuiii_4 : 



N2 ) 



i.parity = 0; 

} 

else 
{ 

if (sil == si2) 
{ 



/* not end of loop */ 



} 

else 
{ 



++inum_l ; 
++inuin_2; 

factor_il = factor_i2 = Sqrt2; 

i_parity =1; /* end of loop && sl=s2 */ 



++inujii_2: 
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++inim_3 ; 

factor_i2 = factor_i3 = Sqrt2; 
i.parity = -1; /* end of loop 



sl!=s2 */ 



element[l] = physix(jm, jt , sj 1 , sj2, im, it , sil , si2) ; 

element[2] = physix(jm, jt , sj 1 , sj2, im, it , -sil , -si2) ; 

element [3] = physix(jm, jt , sj 1 , sj2, im,N2+l-it , si2, sil) ; 

element [4] = physix(jm, jt , sj 1 , sj2, im,N2+l-it , -si2 ,-sil) ; 

element[5] = physix(jm, jt , -sj 1 ,-sj2 , im, it , sil , si2) ; 

element[6] = physix(jm,jt , -sj 1 ,-sj2 , im, it ,-sil ,-si2) ; 

element[7] = physix(jm, jt,-sjl,-sj2,im,N2+l-it,si2,sil) ; 

element[8] = physix(jm, jt,-sjl,-sj2,im,H2+l-it,-si2,-sil) ; 

element[9] = physix(jm,N2+l-jt,sj2,sjl,im,it,sil,si2) ; 
elementdO] = physix(jm,N2+l-jt,sj2,sjl,im,it,-sil,-si2) ; 
element[ll] = physix(jm,N2+l-jt,sj2,sjl, 

im,N2+l-it,si2,sil) ; 
element[12] = physix(jm,N2+l-jt,sj2,sjl, 

im,N2+l-it,-si2,-sil) ; 

element[13] = physix( jm,N2+l-jt , -sj2 , -sj 1 , im, it ,sil , si2) ; 
element[14] = physlx( jm,N2+l-jt , -sj2 , -sj 1 , 

im, it ,-sil , -si2) ; 
element[15] = physlx( jm,N2+l-jt , -sj2 , -sj 1 , 

im,N2+l-it,si2,sil) ; 
element[16] = physix(jm,N2+l-jt , -sj2, -sj 1 , 

im,N2+l-it,-si2,-sil) ; 



/* 



=> take care of phase <= 



*/ 



if (abs(sil+si2-J_z) == 2) 
{ 

element [2] = -element [2] ; 
element [3] = -element [3] ; 
element [14] = -element [14] ; 
element [15] = -element [15] ; 

} 

if (abs(sjl+sj2-J_z) == 2) 
{ 

element [5] = -element [5] ; 
element [8] = -element [8] ; 
element [9] = -element [9] ; 
element [12] = -element [12] ; 



if ((abs(sil+si2-sjl-sj2) ==2) I I 
(abs(sjl+sj2-sjl-sj2) == 3)) 

element [6] = -element [6] ; 
element [7] = -element [7] ; 
element [10] = -element [10] ; 
element [11] = -element [11] ; 



/* Kinetic energy is diagonal */ 

T = 4*(mu[im]*mu[im] + m*m) ; 

if ( (j -parity >= 0) && (i.parity >=0)) 

/* CH = ++ */ 

{ 

mat_l[IndexMat[0]] = ALPHA/4*f actor.j l*f actor.il* 
( 
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element [1] + element [2] - element [3] - element [4] 

+ element [5] + element [6] - element [7] - element [8] 

- element [9] - element [10] + element [11] + element [12] 

- element [13] - element [14] + element [15] + element [16] 
); 

if ( jnum_l == lnuiii_l ) mat_l [IndexMat [0] ] += T; 
++IndexMat [0] ; 



mat_2 [IndexMat [1] ] = ALPHA/4*f actor_j2*f actor_i2* 
/* CH = +- */ 
( 

element [1] - element [2] - element [3] + element [4] 

- element [5] + element [6] + element [7] - element [8] 

- element [9] + element [10] + element [11] - element [12] 
+ element [13] - element [14] - element [15] + element [16] 

); 

if ( jnum_2 == iniim_2 ) mat_2 [IndexMat [1] ] += T; 
++IndexMat [1] ; 

if ( (j.parity <= 0) && (i.parity <=0)) 

/* CH = -+ */ 

{ 

mat_3 [IndexMat [2]] = ALPHA/4*f actor. j3*factor_i3* 
( 

element [1] + element [2] + element [3] + element [4] 
+ element [5] + element [6] + element [7] + element [8] 
+ element [9] + element [10] + element [11] + element [12] 
+ element [13] + element [14] + element [15] + element [16] 
); 

if ( jnum_3 == inum_3 ) mat_3 [IndexMat [2] ] += T; 
++ IndexMat [2] ; 

} 

if ( (j.parity == 0) && (i.parity ==0)) 

/* CH = — */ 

{ 

mat_4 [IndexMat [3]] = ALPHA/4* 
( 

element [1] - element [2] + element [3] - element [4] 

- element [5] + element [6] - element [7] + element [8] 

+ element [9] - element [10] + element [11] - element [12] 

- element [13] + element [14] - element [15] + element [16] 

); 

if ( jnuin_4 == lnum_4 ) mat_4 [IndexMat [3] ] += T; 
++ IndexMat [3] ; 

} 

} /* three column loops: <spin>,<theta>,<mu> */ 

} 

} 

} /* three row loops: <spin>,<theta>,<mu> */ 

} 

} 

time(&time2) ; 

f printf (f _log, " ...matrix: '/,d sec\n" ,time2-timel) ; 

return ; 

} 
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Asymmetrix 

/************************************************************************/ 



/* */ 

/* program : ASYMMETRIX */ 

/* */ 

/* class : subroutine of MESONIX */ 

/* */ 

/* file name : "asymmetrix. c" */ 

/* */ 

/* purpose : ASYMMETRIX generates the Hamilton matrix of the */ 

/* system, but WITHOUT any use of syimnetries . It is */ 

/* as well a test case for the KRAUT code, as a first */ 

/* step to the J_z<>0 case, which has OTHER symmetry */ 

/* properties. It uses the subroutine PHYSIX for */ 

/* calculation of the single elements. */ 

/* */ 

/* */ 

/* */ 
/* External functions: */ 


/* */ 

/* physix subroutine of 'mesonix', calculates one matrix */ 

/* element */ 

/* */ 



void asymmetrixO 

/* create Hamiltonian, no symmetries */ 
{ 

int jm, jt,sjl,sj2; /* loop variables for rows */ 
int im, it , sii , si2; /* loop variables for columns */ 
double T; /* kinetic energy */ 

int j_num,i_num; /* actual Hamiltonian row/column */ 
printfC > ASYMMETRIX\n") ; 



j_num = 0; 
IndexMatrix = 0; 

/* =====================> rows of Hamiltonian <================= */ 



for ( jm=l; jm <= Nl; ++jm) /* loop for <mu> */ 
{ 

for ( jt=l; jt <= N2; ++jt) /* loop for <theta> */ 
{ 

printf ("jm/Nl ="/.2d /"/.2d, jt/N2='/.2d /7,2d\n" , jm.Nl , jt ,K2) ; 

for ( sjl = -1; sjl <= 1; sjl += 2 ) /* loop for <spinl> */ 
{ 

for ( sj2 = -1; sj2 <= 1; sj2 += 2 ) /* loop for <spin2> */ 
{ 

i_num = 0; 
++j_num; 

/* ===============> columns of Hamiltonian <=============== */ 

for ( im=l; im <= Nl; ++im) /* loop for <mu> */ 
{ 

for ( it=l; it <= N2; ++it) /* loop for <theta> */ 
{ 

for ( sil = -1; sil <= 1; sil += 2 ) /* loop for <spinl> */ 
{ 

for ( si2 = -1; si2 <= 1; si2 += 2 ) /* loop for <spin2> */ 
{ 
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++i_num; 

ma[IndexMatrix] = ALPHA*physix( jm, jt , -sj 1 , -sj2, im.it ,-sil,-si2) ; 
if ( j_iium == i_num ) ma [IndexMatr Ix] += 4* (mu [im] *mu [im] + m*m) ; 
++IndexMatrix; /* matrix element counter */ 

} /* four column loops: <spin2>,<spinl>,<theta>,<mu> */ 

} 

} 

} 

} /* four row loops: <spin2>,<spinl>,<theta>,<mu> */ 

} 

} 

} 

return ; 
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Physix 



/* */ 

/* program : P H Y S I X */ 

/* */ 

/* class : sub-subroutine of MESONIX, subroutine of GEHETIX */ 

/* */ 

/* file name : "physix. c" */ 

/* */ 

/* purpose : PHYSIX calculates one single element of the */ 

/* Hamilton matrix of the system. */ 

/* Takes care of singularities with Coulomb trick. */ 

/* */ 
/************************************************************************/ 



double Int(int n) 

/* the formula for the integral over cos(nx)/(a+b cos(x))*/ 
{ 

double vz; 
vz=l . ; 

if (n'/.2==0) vz = -1.0; 

return vz*pow(A, 1 .0-n)*pow(B/k[l] /k[2] ,1.0*n); 

} 

double GKint j, int i, int n) 

/* diagonal matrix element, parallel spins */ 

{ 

double help; 

if ( (k[i]==0) && (k[j]==0) ) help = 0.0; 

else 

{ 

help = m*m*(1.0/x[i]/x[j]+l/(l-x[j])/(l-x[i]))*Int(abs(l-n)) 
+ k [i] *k [ j ] /x [ j] /x [i] / (1-x [ j ] ) / (1-x [i] ) *Int (abs (n) ) ; 

} 

return help; 

} 

double G2(int j, int i, int n) 

/* diagonal matrix element, anti-parallel spins */ 
{ 

double help; 

help = l/x[i]/x[j] + 1.0/(l-x[j])/(l-x[i]); 

help = (m*m*help + k[j]*k[j]/x[j]/(l-x[j] ) 

+ k [i] *k [i] /x [i] / ( 1-x [i] ) ) *Int (abs (n) ) +k [i] *k [ j] * 

(Int (abs (1-n) ) /x [i] /x [ j] +Int (abs ( 1+n) ) / (1-x [i] ) / (1-x [ j] ) ) ; 

return help; 

} 

double G3(int j, int i, int n) 
/* off -diagonal matrix element */ 
{ 

double help; 

if ( k[j]==0 ) help = 0.0; 
else help = -m/x[i]/x[j]* 

(k [i] *Int (abs ( 1-n) ) -k [ j] * (1-x [i] ) / ( 1-x [j ] ) *Int (abs (n) ) ) ; 

return help; 

} 

double G3_star(int j, int i, int n) 
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/* off-diagonal matrix element */ 
{ 

double help; 

if ( k[j]==0 ) help = 0.0; 
else help = m/(l-x[i] )/(l-x[j] )* 

(k [i] *Int (abs (1-n) ) -k [ j ] *x [i] /x [ j ] *Int (abs (n) ) ) ; 

return help; 

} 

double G4(int j, int i, int n) 
/* off-diagonal matrix element */ 
{ 

double help; 

help = -m*m*(x[i]-x[j])*(x[i]-x[j])/(l-x[i])/(l-x[j]) 

/x [i] /x [ j ] * Int (abs (n) ) ; 
return help; 

} 

/* > Functions for dynamical annihilation graph < */ 



double Fl(int j, int 1, int n) 

/* dynamical annihilation graph: parallel spins, diagonal */ 
{ 

double help, omega; 

help=0.0; 

if (abs(n)==l) 

{ 

omega= ( (m*m+k [i] *k [i] ) /x [i] / (1 . 0-x [i] ) 

+ (m*m+k[j]*k[j])/x[j]/(1.0-x[j]))/2.0; 
help=2 . O*m*m/omega* (1 . 0/x [i] +1/ ( 1-x [i] ) ) * ( 1/x [ j ] +1 . 0/ ( 1-x [ j ] ) ) ; 

} 

return help; 

} 

double F2(int j, int i, int n) 

/* dynamical annihilation graph: antiparallel spins, diagonal, I */ 
{ 

double help, omega; 

help =0.0; 

if (n==0) help=4.0; /* seagull graph */ 
if (abs(n)==l) /* dynamic graph */ 

{ 

omega= ( (m*m+k [i] *k [i] ) /x [i] / (1 . 0-x [i] ) 
+ (m*m+k[j]*k[j])/x [j]/(1.0-x [j]))/2.0; 
help = 2.0/omega*k[i]*k[j]/x[i]/x[j] ; 

} 

return help; 

} 

double F2_star(int j, int i, int n) 

/* dynamical annihilation graph: antiparallel spins, diagonal, II */ 
{ 

double help, omega; 

help =0.0; 

if (n==0) help=4.0; /* seagull graph */ 
if (abs(n)==l) /* dynamic graph */ 

{ 

omega= ( (m*m+k [i] *k [i] ) /x [i] / (1 . 0-x [i] ) 
+ (m*m+k[j]*k[j])/x[j]/(1.0-x[j]))/2.0; 
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help = 2.0/omega*k[i]*k[j]/(1.0-x[i])/(1.0-x[j]); 

} 

return help; 

} 

double F3(int j, int i, int n) 

/* dynamical annihilation graph: spin non-diagonal */ 
{ 

double help, omega; 

help=0.0; 

if (abs(n)==l) /* dynamic graph */ 

{ 

omega= ( (m*m+k [i] *k [i] ) /x [i] / ( 1 . 0-x [i] ) 
+ (m*m+k[j]*k[j])/x[j]/(1.0-x[j]))/2.0; 

help = 2.0*m/omega*k[i]/(l-x[i])*(l/x[j]+1.0/(l-x[j])); 

} 

return help; 

} 

double F3_star(int j, int i, int n) 

/* dynamical annihilation graph: spin non-diagonal */ 
{ 

double help, omega; 
help=0.0; 

if (abs(n)==l) /* dynamic graph */ 

{ 

omega=( (m*m+k [i] *k [i] ) /x [i] / ( 1 . 0-x [i] ) 

+ (m*m+k[j]*k[j])/x[j]/(1.0-x [j]))/2.0; 
help = -2.0*m/omega*k[i]/x[i]*(l/x[j]+1.0/(l-x[j])); 

} 

return help; 

} 

double F4(int j, int i, int n) 

/* dynamical annihilation graph: all spins antiparallel, non-diagonal */ 
{ 

double help, omega; 

help = 0.0; 

if (n==0) help=4.0; /* seagull graph */ 
if (abs(n)==l) /* dynamic graph */ 

{ 

omega= ( (m*m+k [i] *k [i] ) /x [i] / ( 1 . 0-x [i] ) + 
(m*m+k[j]*k[j])/x[j]/(1.0-x[j]))/2.0; 
help = -2.0/omega*k[i]*k[j]/x[i]/(l-x[j]); 

} 

return help; 

} 

double annihilation_spintable(int sjl, int sj2, int sil, int si2) 
{ 

double result; 

/* Evaluation of SPIN-TABLE of dynamic annihilation graph */ 
/* J_z == */ 

/* (row-wise from top to bottom, left to right) */ 
/* (final is row index, initial is column index) */ 

result= 0.0; 

if (J_z==-1) /* 'transpose' matrix */ 

{ 

sjl = -sjl; 
sj2 = -sj2; 
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sil = -sil; 
si2 = -si2; 

} 

if (((sjl 1) && (sj2 D) I I ((sil 1) && (si2 1))) return 0.0; 

else 
{ 

if (sjl > 0) /* first and second row */ 

{ 

if (sj2 > 0) /* first row */ 

{ 

if (sil > 0) 
{ 

if (si2 > 0) result = Fl(2,l,J_z); 
else result = F3(2,l,J_z); 

} 

else result = F3_star(2,l, J_z) ; 

} 

else /* second row */ 

{ 

if (sil > 0) 
{ 

if (si2 > 0) result = F3_star(l,2, J_z) ; /* perm. */ 
else result = F2_star(2,l, J_z) ; 

} 

else result = F4(2,l,J_z); 

} 

} 

else /* third row */ 



if (sil > 0) 
{ 

if (si2 > 0) result = F3_star(l,2, J_z) ; /* permuted */ 

else result = F4(l,2,J_z); /* permuted */ 

} 

else result = F2(2,l,J_z); 



} 

return result: 



double general_spintable(int sjl, int sj2, int sil, int si2) 
/* spin table for J_z=n */ 
/* */ 

/* (row-wise from top to bottom, left to right) */ 
/* (final is row index, initial is column index) */ 
{ 

double result ; 

if (sjl > 0) /* first and second row */ 

{ 

if (sj2 > 0) /* first row */ 

{ 

if (sil > 0) 
{ 

if (si2 > 0) result = Gl(2,l,J_z); 
else result = G3_star(2,l, J_z) ; 

} 

else 
{ 

if (si2 >0) result = G3(2,l,J_z); 
else result = 0.0; 

} 

} 

else /* second row */ 
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if (sil > 0) 



else 



if (si2 > 0) result = G3_star (1 ,2 , J_z) ; /* permuted */ 
else result = G2(2,l,J_z); 



if (si2 >0) result = G4(2,l,J_z); 
else result = -03(1 ,2 ,-J_z) ; /* permuted */ 



} 

} 

else /* third and forth row */ 

{ 

if (sj2 > 0) /* third row */ 

{ 

if (sil > 0) 
{ 

if (si2 > 0) result = 03(1,2, J_z) ;/* permuted */ 
else result = 04(2,l,J_z); 

} 

else 
{ 

if (si2 > 0) result = 02(2, 1 , -J.z) ; 

else result = -03_star(l,2,-J_z) ;/* permuted */ 

} 

} 

else /* forth row */ 

{ 

if (sil > 0) 
{ 

if (si2 > 0) result = 0.0; 

else result = -03(2,1,-J_z) ; 

} 

else 
{ 

if (si2 >0) result = -03_star (2, 1 ,-J_z) ; 
else result = 01(2,1,-J_z); 

} 

} 

} 

if (Aimi!=0) result += aiinihilation_spintable(sjl,sj2,sil,si2) ; 
return result; 



int spin2number (int si, int s2) 

/* transforms (uu,ud,du,dd) — > (0,1,2,3) */ 

{ 

return 2-sl-s2-(l-s2)/2; 

} 



double physix(int jm,int jt.int sjl.int sj2,int im.int it, int sil, int si2) 
/* main routine: calculates one Hamiltonian matrix element */ 


/* parameters: (j=final, i=initial) */ 
/* jm.jt = <mu>,<theta> */ 

/* sjl,sj2 = <spin electron>,<spin positron> */ 

{ 

double result, jacobian.epsilon; 

V_seag = -2; /* seagull graph */ 
epsilon = le-8; 
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/* ========> transformation: (mu,theta) — >(x,k) <============ */ 

/* j — > 2 */ 

/* i — > 1 */ 

x[2] = (1.0 + mu[jm] *theta[jt] /sqrt (m*m + inu[jm] *mu [jm] ) ) /2 . ; 
k[2] = mu[jm]*sqrt(1.0-theta[jt]*theta[jt]); 

/* ======================> Coulomb trick <=================== */ 



if ((jm==im) && (jt==it) && (sjl==sil) && (sj2==si2)) 
{ 

result = 0.0; /* default */ 

result = Coulomb_discrete( jm.sjl , sj2) ; /* discrete counterterm */ 

if (old_CT!=0) result += old_continous(jm, jt) ; 

else result += CT[spin2number(sjl,sj2)] [jm-1] [jt-1] ; 

/* continuous counterterm*/ 

} 

/* =====================> No Coulomb trick <================== */ 

else 
{ 

x[l] = (1.0 + mu[im] *theta[it] /sqrt (m*m + mu[im] *mu[im] ) )/2 .0; 

k[l] = mu[ini] *sqrt (1 .0-theta[it] *theta[it] ) ; 



if ( jt==it && jm==iiii ) return 0.0; 

/* Note: aa = -a(Kraut) , 09/11/95 */ 

/* */ 

/* -a,A,B as defined in KPW, etc. */ 

aa = (x[l]-x[2])*(x[l]-x[2])*m*m/2.0*(l/(l-x[2])/(l-x[l]) 
+l/x[l]/x[2]) + k[l]*k[l] + k[2]*k[2] 
+ (x[l]-x[2])/2*(k[l]*k[l]*(l/(l-x[l])-l/x[l]) - 
k [2] *k [2] * (1/ (1-x [2] ) -1/x [2] ) ) ; 

A = 1.0/sqrt( aa*aa - 4*k[l]*k[l]*k[2]*k[2] ); 

B = (1.0-aa*A)/2.0; 



jacobian = sqrt(wmu[jm] *wmu[im] *wtheta[jt] *wtheta[it] ) 

*mu [ jm] *mu [im] *sqrt ( 4*x [1] * (1-x [1] ) *x [2] * ( 1-x [2] ) 
/sqrt(m*m + mu[im] *mu[im] )/sqrt (m*m + mu[jm] *mu[jm] ) ) ; 

result = jacobian/PI*general_spintable(sjl,sj2,sil,si2) ; 

} 

return result; 

} 



Arithmetix 



/* */ 
/* program : ARITHMETIX */ 
/* */ 
/* class : subroutine of MESONIX */ 
/* */ 
/* file name : "arithmetix . c" */ 
/* */ 
/* purpose : ARITHMETIX diagonalizes the four block matrizes */ 
/* of the Hamiltonian and sorts the Eigenvalues from */ 
/* small to big ones. */ 
/* The eigenf unctions are written in a file, if */ 
/* <print_EF>!=0 */ 
/* */ 
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/* */ 

/* External functions: */ 


/* */ 

/* mOlcaf NAGLib routine, sorts vector */ 

/* */ 



extern double m01caf_( double *r, int *ml, int *m2, char *order, 

int *IFAIL); 
/* sorting routine from NAGLib */ 

double normalisationCint i, int j) 

/* the factor to resubstitute Phi — > Psi */ 

{ 

double ai j ; 

aij = mu[i] *mu[i] *(m*m+mu[i] *mu[i] *(1 .0-theta[j] ) ) 

II . 0/ sqrt (m*m+mu [i] *mu [i] ) ; 
aij = 1.0/sqrt(wtheta[j]*wmu[i]*aij) ; 
if (j<(N2+l)/2) aij *= 1.0/sqrt(2.0) ; 
return aij ; 

} 

double asym_normalisation(int i, int j) 

/* the factor to resubstitute Phi — > Psi */ 

{ 

double aij ; 

aij = mu [i] *mu [i] * (m*m+mu[i] *mu[i] *(1 .0-theta[j] ) ) 

/2 . 0/sqrt (iii*m+niu[i] *mu[i] ) ; 
aij = 1 . 0/sqrt (wtheta[j] *wmu[i] *aij ) ; 
aij *= 1. 0/sqrt (2.0) ; 
return aij ; 

} 

void store_asy_eigenfunctions() 

/* write eigenf unctions of asynmetric Hamiltonian in file ' EFunc_asy_<l> . dat ' */ 
{ 

int i,j,l,n; 
char *name_EF [4] ; 
FILE *f p [4] ; 

for (1=0; 1<=3; 1++) name.EF [1] = (char *)malloc (120) ; 

for (n=0; n<=5; n++) /* <n>th eigenvalue */ 

{ 

for (1=0; 1<=3; 1++) 
{ 

sprintf (name.EF [1] , "%sEigenFunctions/EFy,ld_J7,ldasy_7.1d.dat" , 

def ault_directory,n, J_z,l) ; 

printf ("*/.s\n" ,name_EF[l] ) ; 
fp [l]=f open (name.EF [1] ,"w"); 

} 

for (i=l; i<=Nl; i++) 

for (j=l; j<=N2; j++) 
{ 

for (1=0; 1<=3; 1++) 
{ 

x[l] = (1.0 + mu[i]*theta[j]/sqrt(m*m + mu[i] *mu[i] ) )/2. 0; 

k[l] = mu[i]*sqrt(1.0-theta[j]*theta[j]); 

fprintf (fp[l] ,""/.18.12f "/.18.12f '/.IS . 12f \n" ,x [1] ,k[l] , 
V [n*4*Nl*N2+ (i-1) *N2*4+4* ( j -1) +1] * 
asym_normalisation(i, j)) ; 

} 

} 
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for (1=0; 1<=3; f close (fp [1] ) ; 

} /* loop n */ 

} 



void store_eigenf unctions () 

/* write eigenf unctions in files 'EFunc_<n>.dat' */ 
{ 

int i , j ; 

char name.EF [120] ; 
FILE *fpl; 

/* CH = ++ */ 

sprintf (name_EF, "'/.sEigenFunctions/EFunc.l .dat" ,def ault_directory) ; 
fpl=f open(name_EF, "w") ; 
for (i=l; i<=Nl; i++) 

for (j=l; j<=(N2+l)/2; 

{ 

fprintf (fpl,"%18.12f\n",v_l[(i-l)*N2+2*j-2] *normalisation(i, j)) ; 

if (j! = (N2+l)/2) fprintf (fpl,"y.l8.12f\n",v_l[(i-l)*N2+2*j-l]*normalisation(i,j)); 

else fprintf (f pi , "\n") ; /* dummy line */ 

} 

f close(fpl) ; 
/* CH = +- */ 

sprintf (nctme_EF, "y,sEig6nFunctions/EFunc_2 .dat" .default .directory) ; 
fpl=f open(name_EF, "w") ; 
for (i=l; i<=Nl; i++) 

for (j=l; j<=(N2+l)/2; 

{ 

fprintf (f pi , "118 . 12f \n" , v_2 [ (i-1) * (N2+1) +2* j -2] *normalisation(i , j ) ) ; 
fprintf (f pi , "'/.IS . 12f \n" , v_2 [ (i-1) * (H2+l)+2*j -1] *normalisation(i , j ) ) ; 

} 

f close(fpl) ; 
/* CH = -+ */ 

sprintf (name.EF, '"/.sEigenFunctions/EFunc.S.dat" .default .directory) ; 
f pl=f open(name_EF, "w") ; 
for (i=l; i<=Nl; i++) 

for (j=l; j<=(N2+l)/2; j++) 

{ 

fprintf (f pi , "'/.IS . 12f \n" , v_3 [ (i-1) *N2+2* j-2] *normalisation(i , j ) ) ; 

if (j!=(N2+l)/2) 

fprintf (f pi, '"/.IS. 12f\n" , v_3 [(1-1) *N2+2* j-1] *normalisation(i , j ) ) ; 
else fprintf (f pi . "\n") ; /* dummy line */ 
} 

f close(fpl) ; 

/* CH = — */ 

sprintf (name_EF, "/',sEigenFunctions/EFunc_4.dat" ,def ault_directory) ; 
fpl=fopen(name_EF,"w") ; 
for (i=l; i<=Nl; i++) 
{ 

for (j=l; j<(N2+l)/2; j++) 
{ 

fprintf (f pi , "'/.IS . 12f \n" , v_4 [ (i-1) * (N2-1) +2*j -2] *normalisation(i , j ) ) ; 
fprintf (f pi. "'/.IS. 12f\n", v_4[(i-l)*(H2-l)+2*j-l] *normalisation(i, j)) ; 

} 

fprintf (fpl."\n"); /* dummy line */ 
fprintf (fpl."\n"); /* dummy line */ 

} 

f close(fpl) ; 

/* ================> write (x.k) in file 'x_k.dat' <======== */ 
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sprintf (name_EF, "%sEigenFunctions/x_k.dat" ,def ault.directory) ; 
fpl = fopen(name_EF, "w") ; 
for (i=l; i<=Nl; 

for (j=l; j<=(N2+l)/2; 

{ 

x[l] = (1.0 + mu[i]*theta[j]/sqrt(m*m + mu [i] *mu [i] ) ) /2. 0; 

k[l] = mu[i]*sqrt(1.0-theta[j] *theta[j] ) ; 

f printf (fpl , "•/.22 . 16f •/.22 . 16f ,x [1] ,k [1] ) ; 

f printf (f pi , ""/.22 . 16f 7.22. 16f \n" ,x[l] ,k[l]); 

} 

f close (fpl) ; 
/* mu values */ 

sprintf (name_EF, '"/.sEigenFunct ions/abscissae . dat " , default .directory) ; 
fpl = f open(name_EF, "w") ; 

for (i=l; i<=Nl; i++) f printf (fpl , "•/.22. 16f\n" ,mu[i] /p.bohr) ; 
f close (fpl) ; 



void arithmetix(int dim_l, int dim_2, int dim_3, int dim_4) 

/* main function: diagonalize Hamiltonians */ 

{ 

int i; /* loop variable */ 

int IA_nag,N_nag,IV_nag,IFAIL_nag; /* variables for NAGLib */ 

int ml,m2; /* variables for mOlcaf */ 

FILE *ffp; 



time(&timel) ; 

printf ("\n > ARITHMETIX\n") ; 

if ((testa (4096+16384)) == 0) 

/* J_z=0 case, 4 bloc matrices */ 

{ 

/* ==========> diagonalize the 4 bloc matrices <============ */ 

IA_nag = IV_nag = N_nag = dim_l ; IFAIL_nag=-l ; 

f02abf_(mat_l,&IA_nag,&N_nag,r_l,v_l,&IV_nag,e_l, felFAIL.nag) ; 
IA_nag = IV_nag = N_nag = dim_2;IFAIL_nag=-l; 

f02abf_(mat_2,&IA_nag,&N_nag,r_2,v_2,&IV_nag,e_2, felFAIL.nag) ; 
IA_nag = IV_nag = N_nag = dim_3;IFAIL_nag=-l; 

f02abf_(mat_3,&IA_nag,&N_neig,r_3,v_3,ftIV_nag,e_3, felFAIL.nag) ; 

IA_nag = IV_nag = N_nag = dim_4; IFAIL_nag=-l ; 

f02abf _(mat_4,&IA_nag,&N_nag,r_4,v_4,&IV_nag,e_4, aiFAIL_nag) ; 

ffp=fopen( "EVs_l.dat" , "w") ; 
f printf (ffp, "C=+,H=+\n") ; 

for (i=l; i<=dim_l; ++i) f printf (ffp, "'/.IS. 12f\n" ,r_l [i-1] ) ; 
f close (ffp) ; 

ffp=fopen("EVs_2.dat","w") ; 
fprintf (ffp,"C=+,H=-\n") ; 

for (i=l; i<=dim_2; ++i) fprintf (ffp, "'/.IS. 12f\n" ,r_2 [i-1] ) ; 
f close (ffp) ; 

ffp=fopen("EVs_3.dat","w") ; 
fprintf (ffp, "C=-,H=+\n") ; 

for (i=l; i<=dim_3; ++i) fprintf (ffp, "7,18. 12f\n" ,r_3 [i-1] ) ; 
f close (ffp) ; 

f fp=f open("EVs_4.dat" , "w") ; 
fprintf (ffp, "C=-,H=-\n") ; 

for (i=l; i<=dim_4; ++i) fprintf (ffp, "7.18. 12f\n" ,r_4 [i-1] ) ; 

f close (ffp) ; 

for (i=l; i<=dim_l+dim_2+dim_3+dim_4; ++i) 
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/* put all EVs together */ 
{ 

if (i<=dim_l) r[i-l] = r_l[i-l]; 

else 

{ 

if (i<=diiii_l+diiii_2) r[i-l] = r.2 [i-l-dim.l] ; 

else 

{ 

if (i<=diiii_l+diiii_2+diiii_3) 

r[i-l] = r_3[i-l-dim_l-dim_2] ; 
else r[i-l] = r_4[i-l-dim_l-dim_2-dim_3] ; 

} 

} 

} 

if (print.EF > 0) store_eigenfunctions() ; 

} /* end of (test&4096==0) */ 

else 

{ 

if ((test&16384)==0) /* asymmetric Hamiltonian (1 bloc) */ 
{ 

printfC (asymmetric case)\n"); 

IA_nag = IV_nag = N.nag = dim_l; IFAIL_nag = -1; 

f 02abf _ (ma , &I A_nag , &N_nag , r , v , &I V_nag , e , &IFAIL_nag) ; 
if (print_EF > 0) store_asy_eigenf unctions () ; 

} 

else /* use C-symmetric Hamiltonian (2 blocs) */ 
{ 

printfC (symmetric case)\n"); 

IA_nag = IV_nag = N_nag = dim_l; IFAIL_nag = -1; 

f 02abf _ (mat_l , &IA_nag, &N_nag , r_l , v_l , felV.nag , e_l , ftlFAIL.nag) ; 

IA_nag = IV_nag = N_nag = dim_2; IFAIL_nag = -1; 

f 02abf _ (mat_2 , felA.nag, &N_nag , r_2 , v_2 , felV.nag , e_2 , &IFAIL_nag) ; 

/* ============> write all eigenvalues in <r[i]> <========== */ 

for (i=0; i<=dim_l-l; ++i) 
{ 

r[i] = r_l[i]; 

printf ("H_EW['/.2d] = ■/.12 .8f \n" , i ,r_l [i-1] ) ; 

} 

for (i=0; i<=dim_2-l; ++i) 
{ 

r[dim_l+i+l] = r_2[i]; 

printf ("H_EW[y.2d] = •/.12.8f \n" , i ,r_2 [i-1] ) ; 

} 

} 

} 

order = "A"; ml = 1; m2 = dim_ l+dim_2+dim_3+dim_4 ; 

mOlcaf _(r,feml ,&m2, order ,&IFAIL_nag) ; /* sort eigenvalues */ 

time(&time2) ; 

fprintf (f_log," . . .diagonalization: °/,d sec\n" ,time2-timel) ; 

return ; 
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Appendix H. Listing of the computer code 



Publix 



/* 








*/ 


/* 


program : 


PUBLIX 




*/ 


/* 








*/ 


/* 


class : 


subroutine of 


MESONIX 


*/ 


/* 








*/ 


/* 


file name : 


"publix. c" 




*/ 


/* 








*/ 


/* 


purpose : 


PUBLIX is the 


output routine of MESOHIX. 


*/ 


/* 








*/ 



void publixO 

/* write eigenvalues in file 'EVs_J<n> . dat ' */ 
{ 

int i; /* loop variable */ 

double bs.bt.chf; /* binding coefficients */ 

FILE *fev; /* pointer onto file */ 

printfC > PUBLIX\n"); 

if ((test & 4) != 0) /* write EVs of 4 blocs in different files */ 
{ 

printf("»»> INFO : 2nd test bit set ON => eigenvalues in 'EV_<n>.dat'\n") ; 

fev = fopen("EV_l.dat","w"); 

fprintfCfev, "Sector: CH = ++\n \n\n" , i ,r_l [i] ) ; 

for (i=0; i<=Nl*N2-l; ++i) fprintf (fev, "'/.IS. 12f \n" ,r_l [i] ) ; 
f close (fev) ; 

fev = f open("EV_2.dat" , "w") ; 

fprintf (fev, "Sector: CH = +-\n \n\n" ,i,r_l [i] ) ; 

for (i=0; i<=Nl*N2+Nl-l ; ++i) fprintf (fev, "7,18. 12f\n" ,r_2 [i] ) ; 

f close (fev) ; 

fev = fopen("EV_3.dat","w"); 

fprintf (fev, "Sector: CH = -+\n \n\n" ,i,r_l [i] ) ; 

for (i=0; i<=Nl*N2-l; ++i) fprintf (fev, ""/.18. 12f\n" ,r_3 [i] ) ; 
f close (fev) ; 

fev = fopen("EV_4.dat","w"); 

fprintf (fev, "Sector: CH = — \n \n\n" , i ,r_l [i] ) ; 

for (i=0; i<=Nl*N2-Nl-l ; ++i) fprintf (fev, "'/.IB. 12f\n" ,r_4 [i] ) ; 
f close (fev) ; 

} 

fev = f open(name_EV, "a") ; /* open file for eigenvalues */ 
i=-l; 

while ( (++i<4*Nl*N2)&&(i<nuinber_EV) ) fprintf (f ev, "'/.IS . 12f \n" ,r [i] ) ; 
for (i=4*Nl*N2+l; i<=number_EV; ++i) fprintf (fev, "\n") ; 
/* fill with empty lines up to <number_EV> */ 

f close (fev) ; 

/* ==============> calculate binding coefficients <============== */ 

bs = 4.0/ALPHA/ALPHA*(2.0-sqrt(r[0] )) ; 

bt = 4.0/ALPHA/ALPHA*(2.0-sqrt(r[l])) ; 

chf = (sqrt(r[l])-sqrt(r[0])) /ALPHA/ALPHA/ALPHA/ALPHA; 

fprintf (f .log," \n B_s = ■/.12.8f \n" ,bs) ; 

fprintf (f .log," B.t = •/.12.8f \n" ,bt) ; 

fprintf (f. log," C.hf = •/.12.8f\n" ,chf ) ; 

return; 

} 
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